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PREFACE 


This technical report presents the development of a three- 
dimensional finite element fracture mechanics analysis of composite 
materials, sponsored in part by NASA-Langley Research Center under 
Research Grant NAG-1-454. This grant was initiated in March 1984, 
preliminary development work having been underway as an in-house effort 
at the University of Wyoming for some time prior to that date. The 
NASA-Langley Technical Monitor was Dr. John H. Crews, Jr. 

All work was performed by the Composite Materials Research Group 
(CMRG) within the Department of Mechanical Engineering at the University 
of Wyoming. The Principal Investigator was Dr. Donald F. Adams, 
Professor of Mechanical Engineering. He was assisted by Mr. Jayant M. 
Mahishi, Ph.D. graduate student in Mechanical Engineering. 

Use of names of manufacturers in this report does not constitute 
official endorsement of such manufacturers, either expressed or implied, 
by the National Aeronautics and Space Administration. 
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SECTION 1 


INTRODUCTION 

Much of the further development and application of fiber- 
reinforced composites in high performance structures now depends mainly 
on the fracture characterization of these materials. Unfortunately, the 
fracture behavior of these materials differs considerably from that of 
homogeneous isotropic materials. 

Even though the mechanics of initiation and growth of cracks in 
composites is not well understood, it is generally believed that micro- 
flaws in the form, for example, of broken fibers, fiber-matrix interface 
debonds, and matrix cracks that exist in composites for various reasons, 
influence the fracture process as much as material inhomogeneity and 
anisotropy. Energy absorption during crack propagation in a multi- 
layered composite laminate is typically due to intralaminar transverse 
cracking, interlaminar delamination, fiber breaks, matrix cracks, 
fiber-matrix debonds, fiber pullouts and matrix yielding (see Figure 1) . 
These failure modes depend on a number of factors such as fiber orienta- 
tion, individual layer thickness, and the constitutive relations that 
describe the mechanical properties of the fiber, the matrix and the 
interface. 

The large number of variables involved in the geometry of com- 
posites makes it very difficult to completely characterize their frac- 
ture behavior experimentally. An analytical model which can represent 
the various aspects of the intrinsic physical failure process is thus 
highly desirable. The main obstacle to deriving such a single analy- 
tical model is one of scale. The fiber breaks, matrix cracks and 
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Figure 1. Energy Absorption Mechanisms During Crack Propagation in a 
Multilayered Fiber-Reinforced Composite Laminate. 




fiber-matrix debonds are of the order of one fiber diameter in size, 
while the intralaminar transverse cracks and interlaminar delaminations 
are several orders of magnitude larger than the fiber diameter. Accord- 
ingly, two different analytical approaches are commonly used, viz., the 
macromechanics approach and the micromechanics approach. The macro- 
mechanics approach, which treats the composite as a homogeneous aniso- 
tropic continuum, has been relatively successful in predicting some 
macro scale failure modes. However, it does not account for failure 
modes associated with the material heterogeneity. On the other hand, the 
micromechanics approach has been very useful in studying the influence 
of microflaws on the fracture behavior of composites. Obviously a 
satisfactory approach is one that combines both the micromechanics and 
the macromechanics approaches . 

The main objective of the present work was to develop a fracture 
criterion for composites based upon a combined micromechanics/macro— 
mechanics analysis. The historical development and present state of the 
art of fracture mechanics of composites are reviewed in Section 2. The 
theory associated with the proposed integrated micromechanics and 
macromechanics fracture criterion (IMMFC) for fiber-reinforced com- 
posites is explained in Section 3. A critical elastic strain energy 
release rate in the presence of plasticity has been defined and used as 
a criterion for crack initiation and propagation in both the micro- 
mechanics and the macromechanics analyses. The energy release rate 
criterion and the finite element technique of crack growth simulation in 
three-dimensional finite element models based on the virtual crack 
extension method are explained in Section 4. Preliminary results using 
the micromechanical analysis are given in Section 5. Graphite/epoxy 
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unidirectional composite models are used to study the initiation and 
propagation of microcracks from microscopic flaws. The fracture tough- 
ness values, which represent the energy release rate at the onset of 
fast fracture under different loading conditions, are evaluated using 
the crack growth resistance curve method. These values are subsequently 
used as critical energy release rates in the integrated analysis. 

In Section 6, example results using the macromechanical analysis 
are presented. The initiation and propagation of a delamination crack 
in central notched and single-edge notched [±45/0] g graphite/epoxy 
laminates subjected to inplane tensile stresses normal to the notch are 
studied. The delamination at ±45 and -45/0 interfaces, and material 
yielding at different stages of loading, are presented. 

Section 7 contains the application of the integrated micro- 
mechanical and macromechanical fracture criterion to predict the onset 
and growth of cracks in a single-edge notched [±45/0] graphite/epoxy 
composite laminate. A qualitative comparison with experimental results 
is also presented. 

Section 8 contains a summary, conclusions and suggestions for 
future work. 
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SECTION 2 


BACKGROUND REVIEW 
2.1 Classical Fracture Mechanics 

Fracture mechanics is primarily concerned with the strength of 
materials in the presence of cracks. Until the early twentieth century, 
a satisfactory explanation had not been found for experimental obser- 
vations that fracture occurs at a much lower stress level than that 
needed to separate adjacent atoms across the fracture surface. Griffith 
[1,2] proposed a theory based upon an energy balance associated with the 
fracture process, equating the release of elastic potential energy 
during crack extension under constant load point displacement to the 
surface energy of the newly created crack surface. Realizing that a 
certain amount of plastic deformation takes place near the crack tip in 
most engineering materials, Irwin [3] and Orowan [4] modified the 
Griffith energy balance equation by equating the energy release rate to 
the plastic energy dissipation rate and the surface energy absorption 
rate, which formed the basis for linear elastic fracture mechanics 
(LEFM) . 

Irwin [5] postulated that the critical energy release rate or 
"fracture toughness" can be related to the stress intensity factor, 
which is a measure of the strength of the elastic singularity of the 
crack tip strain field. Such a fracture characterization based on LEFM 
theory is applicable only to those materials in which plastic yielding 
is confined to a very small region around the crack tip, i.e., when the 
stress field associated with the crack tip is mainly elastic. 



In the last two decades, a considerable amount of work has been 
done in developing methods to characterize fracture associated with 
medium to large scale plastic yielding in ductile materials. The 
plastic zone correction method by Irwin, et. al. [6], the J-integral 
method by Rice [7] and Cherepanov [8], the crack opening displacement 
(COD) method by Wells [9] and Cottrell [10], the near tip strain criter- 
ion by Ke and Liu [11], and the equivalent energy method by Witt [12] 
are some of the important theories that have been proposed to deal with 
crack tip plasticity in the absence of crack growth. The crack growth 
resistance curve (R-Curve) method [13], the nonlinear energy method by 
Liebowitz and Eftis [14], the strain energy density criterion by Sih 
[15], and the crack tip opening angle (CTOA) criterion by de Koning [16] 
and Shih, et. al. [17] are some of the theories proposed for ductile 
fracture associated with subcritical crack growth. 

Irwin [6], in an early attempt, accounted for plastic yielding near 
the crack tip in the classical elastic analysis by assuming a slightly 
larger than actual crack size. This plastic zone correction method 
proposed by Irwin is essentially an empirical correction to the linear 
elastic solution. 

Rice [7] and Cherepanov [8] independently introduced a path- 
independent integral (referred to as the J-integral) , which is expressed 
as 

J = / (wdy - a. .n.u. )ds (1) 

J i] J i.x 

r 

where w is the strain energy density, r is any curvilinear path sur- 
rounding the crack tip, starting from any point on the lower crack 
surface and ending on the upper crack surface, traversing in a 
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counterclockwise direction, n^ is the outward unit normal vector to the 

curve, a., is the traction vector parallel with n., u. is the 
ij j i,x 

derivative of the displacement vector, and s is the arc length. It has 

been shown [18,19] that the J- integral represents the energy release 

rate, and that its value at the onset of cracking (J ) can be used as a 

c 

fracture criterion. 

The concept of crack opening displacement (COD) introduced by Wells 
[9] and Cottrel [10] was successful in taking into account the consider- 
able amount of plastic flow associated with fracture in high toughness 
materials. For the fracture of elastic-perf ectly plastic materials in 
plane stress, the can be energy release rate G related to the crack 
opening displacement 6 as 


G = o y 6 (2) 

where a y is the yield stress of the material in simple tension. Hence 
the COD can be used as an effective fracture criterion. In recent 
studies of elastic-plastic fracture mechanics, the COD has been found to 
be a valuable concept to describe stable crack growth. 

The near tip strain criterion proposed by Ke and Liu [11] is based 
on the fact that the transverse strains near a crack tip at the onset of 
fracture are independent of the specimen geometry and type of loading. 
It has been argued [20] that the plastic strains at all points closer 
than a certain small characteristic distance from the crack tip attain a 
value equal to or greater than a critical value; this distance can be 
treated as a fracture criterion for growing cracks. Experimentally it 
is difficult to measure the strains at the midplane of a plate. However, 


7 



available experimental results based on surface strains [ 11 ] indicate 
significant promise in using strain as a fracture criterion for plates 
of all thicknesses. 

The crack growth resistance curve (R-Curve) concept [13], in which 
the resistance to crack growth is plotted against crack extension, has 
been successfully used to explain the slow and stable crack growth in 
ductile materials which precedes unstable crack growth. The R-Curve, 
which has been shown to be independent of the initial crack size [13], 
represents the energy balance equation G = R, where G is the energy 
release rate and R is the energy absorption rate. The point of crack 
instability occurs when 

dG _ dR ( 3 ) 
da da 

The nonlinear energy method introduced by Liebowitz and Eftis [14] 
takes into account the nonlinearity of the compliance curve of a non- 
linear elastic material. A new fracture parameter G, which is con- 
sidered to be more accurate than the linear elastic fracture toughness 

G , is defined as 
C 


G = 


CG, 


(4) 


where C is a constant determined by the geometry of the nonlinear 
compliance curve. The parameter G is applicable to subcritical crack 
growth as well. It has been claimed that the G value corresponding to 
the onset of stable fracture is independent of the geometry effect of 
the material (including the thickness), whereas G c corresponding to the 
onset of rapid fracture depends on the thickness of the material. 
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The strain energy density theory originally proposed by Sih [15] 

for linear elastic fracture was later extended to explain ductile 

fracture and slow stable crack growth as well [21], According to this 

theory, the crack extension is postulated to occur in the direction of 

minimum strain energy density when the strain energy density factor S 

reaches a critical value S . 

C 

The strain energy density factor S is given by 


S = 


/du. 

r c<d7> 


(5) 


where (— ) is the strain energy density and r is the radius of the core 

region surrounding the crack tip, which is a characteristic of the 
material . 

Shih [17] and de Koning [16] have suggested that the crack tip 
opening angle (CTOA) , measured as the angle subtended by the opened 
surface at a small but fixed distance behind the tip of the growing 
crack, is a satisfactory criterion for stable crack growth. Several 
investigators [17,20] have produced supporting theoretical and numerical 
results. 

The theories discussed above identify certain parameters as a 
criterion for the onset and growth of cracks in ductile materials. 
However, each one has its own merits and demerits relative to the 
others. For materials exhibiting extensive plastic deformation, the 
obvious difference from elastic response is the irreversible strain 
followed by linear elastic unloading. The later feature of ductile 

fracture was not given much importance in the development of elastic- 
plastic fracture theories until recently, when Turner [22,23] noted that 
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the energy release rate in the presence of limit load plasticity is not 
G, but a term he called I, which is greater than G. The parameter I can 
be used as a criterion for the initiation and slow stable crack growth 
in ductile materials. 

2.2 Fracture Mechanics of Composite Materials 

The failure modes associated with fracture in fiber-reinforced 
composites differ considerably from those of homogeneous isotropic 
materials. The microscopic flaws pre-existing in composites for various 
reasons, typically in the form of broken fibers, matrix cracks and 
debonded fiber-matrix interfaces, influence the fracture process as much 
as material heterogeneity and anisotropy. The failure modes associated 
with fracture in composites are typically in the form of transverse 
cracking, delamination, fiber breaks, matrix yielding, matrix cracking, 
and fiber-matrix debonding (as previously indicated in Figure 1). As a 
result, the crack growth often does not occur in a self-similar fashion 
as in the case of isotropic materials. Obviously, the theories of 
classical fracture mechanics cannot be directly applied to composite 
materials. 

Owing to the importance and complexity of the subject, a very large 
volume of literature dealing with different aspects of fracture in 
composites has accumulated within the past two decades. A number of 
theories have been proposed. The various approaches adopted for fracture 
characterization of composites can be broadly classified as 

i) micromechanics approaches 

ii) macromechanics approaches 

iii) anisotropic continuum approaches 
These approaches are discussed in more detail in the following sections. 
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2.2.1 Micromechanics Approaches 


These approaches were primarily developed to estimate the mechan- 
ical properties of unidirectional composites based on the properties of 
their constituent materials. The early development of simple self- 
consistent mechanics models [24-29], variational models [30-34] and 
exact models [35-44] are discussed in References [45-47]. The "material 
model" approach developed by Hedgepath and Van Dyke [48] for the stress 
concentration around a single broken fiber, and the associated matrix 
damage in a unidirectional composite, has been extended by Goree and 
Gross for the case of an arbitrary number of broken fibers [49], as well 
as for longitudinal yielding and splitting of the matrix [50], in a 
three-dimensional model containing equally spaced unidirectional fibers. 
There are a number of papers dealing with cracks at interfaces between 
different materials. Some of these directly address the situation in 
composites [51,52]. In Reference [51], the stress intensity factors for 
a crack running into and crossing an interface are evaluated. The 
elastostatic interaction between a penny-shaped crack and elastic fibers 
has been discussed in Reference [52]. The effects of cracks and imper- 
fections of the fiber-matrix interface contact surface on composite 
properties and the onset of brittle fracture are addressed in Reference 
[53]. 

In the fracture characterization of composites, even more important 
than a rigorous stress analysis [48-53] is the estimation of the energy 
absorption associated with the failure at microflaws. Simple estimations 
of the energy absorption during fiber pullout, fiber-matrix debonding, 
"stress relaxation" due to fiber breaks, "crack bridging" (a process in 
which the fibers are left intact as the crack propagates), and plastic 
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deformation of the matrix have been compiled by Phillips and Tetelman 
[54] and Cooper [55]. The different failure mechanisms and the assoc- 
iated energy absorption calculations have been reviewed in References 
[56,57]. 

By introducing a crack propagation capability into an elasto- 
plastic, generalized plane strain finite element micromechanics . model, 
the energy absorption mechanisms during crack propagation in a region of 
a broken fiber have been studied in References [58,59]. A "failed 11 
element approach developed by Adams [60-64] for crack propagation in 
two-dimensional finite element models was used in these analyses. 

The initiation and propagation of microcracks from broken fibers 
and matrix cracks in a single-fiber model composite has been studied by 
Mahishi and Adams [65-67] using an axisymxnetric, elastoplastic finite 
element analysis. Similar single-fiber models have been used by several 
other authors, e.g., by Ko [68] to study stress concentrations due to 
the broken fiber, and by Gradin, et. al. [69] to study the debonding 
between the fiber and the matrix. Mahishi and Adams observed in their 
studies [65-67] that under fixed grip conditions, the microcracks 
originating from a broken fiber end or from a matrix crack grow in a 
slow and stable manner initially to a certain distance, before growing 
rapidly in an unstable manner across the model. The situation of slow 
and stable subcritical crack growth being similar to that observed in 
the case of crack growth in thin metallic plates led to their use of 
crack growth resistance curves (R-curves) to predict the point of crack 
instability in composites. The energy release rate corresponding to the 
point of instability has been defined as a measure of fracture toughness 
of the model composite. This is a very important development in the 
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micromechanics approach to fracture characterization of composites as, 
in addition to being able to predict a single fracture parameter 
(fracture toughness), it can be used to compare the severity of differ- 
ent types and sizes of microflaws. The axisymmetric model also has the 
advantage of being useful in correlations with experimental results. The 
requisite experiments are simple and easily performed. Such experiments 
are underway at a number of research laboratories [70], including within 
the Composite Materials Research Group at the University of Wyoming. 
Mahishi and Adams [67] have also used the axisymmetric model analysis to 
study the influence of weak fiber-matrix interfaces, curing-induced 
residual thermal stresses, and environmental moisture absorption on the 
fracture behavior of graphite/epoxy model composites. A three-dimen- 
sional finite element model has also been developed by these authors to 
study the stress state near a broken fiber in a unidirectional boron/ 
aluminum composite [71]. The fibers were assumed to be packed in a 
regular square array; a repeating cell was identified consisting of a 
single broken fiber surrounded by an array of continuous fibers. This 
model has considerable flexibility in representing different types of 
microflaws and loading conditions. 

By means of such rigorous three-dimensional finite element analyses 
it is possible to study the effects of sizes and density of different 
types of microflaws on the fracture toughness of a unidirectional 
composite. It is important to verify the predicted results of the 
analysis with experimental data. However, the reduction of experimental 
fracture mechanics data to identify different forms of microcracks is 
not straightforward. Beaumont and Anstice [72] have presented a statis- 
tical approach for the failure analysis of micromechanical fractures in 
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graphite fiber and glass fiber composites. Using a simple micromechanics 
model, they were able to estimate the energy dissipated during partial 
debonding of the fiber-matrix interface, during fiber fracture and 
during fiber pullout. 

A two-dimensional, micromechanical finite element analysis has been 
used by Mandel, et. al. [73] to study the initiation of cracks in a 
steel-fiber-reinforced methacrylate polymer matrix, the results being 
compared with experimental data. Good agreement was obtained between the 
predicted and experimental values of loading required for the initiation 
of cracks in the matrix and subsequent crack arrest by the fibers. 
Badar, et. al. [74] have studied the micromechanisms of fracture in 
short— fiber— reinforced thermoplastics. Williams and Reifsnider [75] have 
used three-dimensional micromechanics models to evaluate the internal 
stress field in applying their strain energy release rate approach to 
the prediction of failure modes in composites. Wells and Beaumont [76] 
have used simple micromechanics models to calculate fracture energy 
associated with fiber debond and fiber pullout mechanisms. A general 
survey is presented by Hashin [77]. 

2.2.2 Macromechanics Approaches 

A distinction should be made between different analytical ap- 
proaches dealing with macroscale cracks in composite laminates, depend- 
ing upon whether the macrocracks considered are at the lamina (ply) 
level or at the laminate level. The interlaminar delamination and 
intralaminar transverse cracking which occur at the lamina level are 
usually treated by idealizing the individual laminae as homogeneous 
anisotropic continua. The physical properties of the individual ideal- 
ized laminae are thus represented by average values, based upon fiber 
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and matrix properties. In the present discussion, this approach to 
macroscale cracks in composites is classified as a macromechanics 
approach. 

At the laminate level a phenomenological fracture analysis has been 
developed to study through- the-thickness cracks, by assuming the entire 
laminate as a homogeneous anisotropic continuum. This approach is 
discussed in the next section. 

The delamination mode of crack growth, which is characteristic of 
multilayered composite laminates, has been attributed to the existence 
of interlaminar stresses near the free edges of the laminate [78,79]. A 
linear elastic stress analysis of the free edge in a composite laminate 
[80] suggests that the interlaminar stresses become singular near the 
free edges, and the signs of these stresses depend on the stacking 
sequence of the laminate [79-81]. Hence some laminates may be more prone 
to delamination than others depending on how the plies are stacked 
together. A review of interlaminar stress effects is given in Reference 
[82]. 

The singular behavior of the interlaminar stresses precludes use of 
any failure criterion based on maximum stress for predicting the onset 
and growth of the delamination. On the other hand, the homogeneous 
continuum idealization enables one to use classical fracture mechanics 
concepts. The critical energy release rate criterion has been used by 
Rybicki, et. al. [83] to predict the stable delamination crack growth in 
composite laminates. They used a finite element analysis and a numerical 
technique for evaluating energy release rates based on the Irwin crack- 
closure integral [84], a technique developed earlier by Rybicki and 
Kanninen [85]. This numerical technique has also been used by several 
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other investigators to characterize delamination crack growth in com- 
posite laminates under static tension, static compression, and fatigue 
loading conditions [86-103]. 

Wang and Crossman [86,95] have further developed the technique, to 
include transverse cracking as well. The principal assumptions in their 
energy release rate approach [86] are that the edge delamination in- 
volves only matrix-dominated fracture, which is assumed to be elastic 
and brittle, and that the crack surface is parallel to the ply inter- 
face. Even though the extent of delamination predicted by this energy 
release rate approach is in agreement with experimental observations 
[87-89, 97-100], there exists a difference between the experimentally 
measured energy release rate values and the predicted values. This is 
possibly because in reality the delamination crack does not propagate in 
a self-similar fashion parallel to the ply interface, but rather takes a 
zig-zag path. There is also some concern about possible material yield- 
ing in the vicinity of the crack tip. 

Delamination crack growth in composite laminates is essentially a 
three-dimensional problem. A quasi three-dimensional finite element 
analysis was used in References [93,99,100]. A fully three-dimensional 
finite element analysis, in which the delamination growth surface can be 
fully simulated, has been used in Reference [91]. Wang and his assoc- 
iates [104-107] have used singular finite element, hybrid-stress models 
to study delamination crack growth in composites under both static and 
cyclic loadings. They used a mixed-mode failure criterion for crack 
growth. The energy release rate was calculated directly from the rela- 
tion between stress intensity factors and energy release rates. 
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The transverse cracking mechanism in composite laminates has also 
been extensively studied [86,95]. Transverse cracks are caused by the 
in-plane tensile stress normal to the fiber direction in a unidirec- 
tional ply. The irregular spacing intervals of the multiple transverse 
cracks observed in experiments are attributed to the presence of micro- 
flaws. The energy release rate approach developed for delamination crack 
growth has been used to study transverse crack growth also [86,95]. 

While considerable progress has been made in understanding the 
delamination and transverse cracking mechanisms in composites incorpor- 
ating brittle or quasi-brittle matrix materials, the problem of tough- 
ened polymer matrix or metal matrix composites, in which large scale 
yielding is associated with the cracking, requires special attention. 
2.2.3 Anisotropic Continuum Approach 

This is a phenomenological approach in which the composite is 
assumed to be a homogeneous anisotropic continuum. Perhaps the success 
of the homogeneous anisotropic continuum idealization in the stress and 
deformation analysis of composites led to the development of this 
approach for the fracture characterization of composites. But the fact 
is that the composite stresses and deformations are averaged properties, 
whereas failure is a localized process in which the material hetero- 
geneity plays an important role. The single most important justification 
for such an idealization may be that the classical fracture mechanics 
concepts can be readily employed. 

The anisotropic elasticity solutions [108,109] have been used in 
References [110-112] to evaluate the stress intensity factors for cracks 
of different sizes and orientations under various loading conditions. 
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A number of fracture mechanics theories have been proposed based on 
the LEFM concept [113,114]. These theories have been reviewed exten- 
sively in References [56,115]. Some of the more important of these 
theories are briefly discussed here. 

Waddoups, et. al. [116] adapted the Bowie [117] solution for a 
crack emanating from a hole in an isotropic material, and derived an 
empirical expression for a crack of length 2a in an infinite body under 
tensile load normal to the crack, i.e., 

1/2 

K x = airOfc + a) x/ * (6) 


where l is the characteristic length of an "intense energy region" at 
the crack tip. The critical stress for crack extension is then 


K 


a c = 


IC 


tt(£ + a) 


1/2 


(7) 


The two parameters K and £ were evaluated from experimental data. The 

XU 

predicted strengths were in good agreement with experimentally measured 
values. Cruse [118] represented a hole of radius R as a crack of half 
length a. The value of a was obtained through comparison of the hole and 
crack solutions for an orthotropic material. The strength of the lami- 
nate was then obtained using LEFM and an experimentally determined value 
of Whitney and Nuismer [114], using a stress concentration ap- 

proach, have proposed two stress criteria for the strength of notched 
composites. The first criterion assumes that fracture occurs when the 
normal stress perpendicular to and ahead of the crack reaches the 
strength of the unnotched laminate at a specific distance d 0 from the 
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crack tip. The distance d 0 is assumed to be a material property, inde- 
pendent of the laminate geometry and stress distribution. In the second 
criterion it is assumed that failure occurs when the average stress over 
a distance a 0 ahead of the crack tip reaches the unnotched tensile 
strength of the laminate. Reasonably good agreement was obtained with 
the limited experimental results available. 

Pipes, et. al. [119], generalizing the point stress and average 
stress criteria of Whitney and Nuismer [114], have introduced a three- 
parameter fracture criterion. A fictitious crack model has been pre- 
sented by Backlund [120]. 

All of the theories discussed so far assume self-similar crack 
growth. In reality, cracks in composite laminates seldom grow in a 
self-similar fashion. This factor alone makes the application of these 
theories for the fracture characterization of composites somewhat less 
effective . 

Harrison [121], to remove the restriction of self-similar crack 
growth, postulated different energy release rates for crack growth in 
the plane of the crack (denoted as G^) and for crack growth normal to 
the crack (denoted as G^) . In his study of splitting in fiber- 
reinforced materials, Harrison gave the condition of splitting as 


G R 

-* < 

G R 

y y 


( 8 ) 


where and are the critical energy absorption rates for crack 
growth in the two directions. More general theories have been developed 
by Sih and Chen [122] and Wu [123,124]. The method of Sih and Chen is 
based on a strain energy density fracture criterion, which not only 
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predicts fracture but also the direction of fracture. This theory also 
introduces a parameter r , the radius of the core region in Eq. (5), 
which can be evaluated analytically or experimentally. The technique has 
been applied both to unidirectional composites and angle-ply laminates. 
Lukshminarayana, et. al. [125] have applied the same technique to cracks 
emanating from circular holes in cross-ply laminates. Experimental 
correlations were limited to transverse cracks running parallel to the 
fibers, but were in good agreement. 

Wu’s theory [123,124] also predicts fracture and the direction of 

fracture. The theory involves application of a failure criterion to the 

region ahead of a macroscopic crack in an anisotropic body. Wu assumes 

that the composite is a homogeneous anisotropic continuum containing 

randomly distributed microscopic cracks. For such configurations he 

further assumes that there exists a small but finite volume of dimension 

r that fully encapsulates a microscopic crack, such that the singular 
c 

stresses are contained within the critical volume and the stresses 
external to the critical volume are bounded. He then postulates that 
failure occurs when the stress vector (o acting on the outer surface of 
the critical volume reaches the value of the strength vector defining 
the failure surface for the material failure. Wu obtained excellent 
agreement with experiments on a Scotchply 1002 glass/epoxy laminate for 
mixed mode loading with cracks parallel to the fibers. In this theory, 
the failure surface must be obtained from experimental studies to 
determine remote properties. 

There are other approaches which deal only with the specific mode 
of crack growth. One example is that of Kulkarni and Rosen [126], which 
is limited to crack growth normal to the crack plane, as observed for 


20 



both unidirectional and general laminates that contain 0° plies. This 
model was originally developed by Zweben [127], who used the detailed 
study of plasticity and crack blunting, in an orthotropic body given by 
Tirosh [128]. In their model, Kulkarni and Rosen [126] have taken the 
region adjacent to the notch as a shear stress transfer region. The 
regions adjacent to the crack are subdivided into a region of shear 
stress transfer, a region of stress concentration, and a region of shear 
stress transfer in the average material. In the application of this 
approach to fracture normal to the original crack surface, it is neces- 
sary to determine two arbitrary parameters experimentally which define 
the size of the region of intense energy adjacent to the flaw, and the 
demarcation between Mode I and Mode II behavior. 

In summary, the different approaches attempting to extend the 
classical fracture mechanics concepts to composite materials are, as 
noted by Kanninen, et. al. [115], nothing but curve-fitting techniques 
with parameters which can be adjusted to fit any experimental data. It 
can be concluded that the fracture behavior of composites is far more 
complex than can be modeled by LEFM. The main reason for the complexity 
is the material heterogeneity. A more general approach to study the 
fracture mechanics of composites should include the affects of material 
heterogeneity. 

Kanninen, et. al. [115,129] have developed an approach which 
combines a ’ micromechanical failure analysis with a macromechanical 
continuum analysis. At the crack tip, a local region is defined consist- 
ing of fibers, the surrounding matrix, and the associated interfaces. 
This local heterogeneous region (LHR) is considered to be surrounded by 
a homogeneous orthotropic continuum. The three constituent materials 
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(which include an interface material) are modeled using "spring-like" 
elements. The size of the LHR is assumed to be large enough such that 
its boundary displacements are given by continuum analyses, and small 
enough relative to the crack length that the singular behavior at the 
crack tip dominates. By varying the fracture properties, Kanninen, et. 
al. [129] were able to demonstrate the occurrence of axial splitting, 
matrix crazing, matrix bridging and fiber bridging. There was also 
qualitative agreement with the experimental results of Brinson and Yeow 
[130] on edge notched unidirectional graphite/epoxy laminates. Represen- 
tation of the local heterogeneous region in a multilayered laminate, and 
in the case of continuous crack propagation, presents some practical 
limitations. 

From this literature study of different approaches to fracture of 
composites it is clear that the micromechanics approach, which deals 
with microscopic failure processes, and the macro mechanics approach, 
which deals with macroscale delaminations and transverse cracking 
failure processes, are very well developed. A conglomeration of both 
microscopic and macroscale failure processes represent the actual damage 
in a multilayered composite. Thus, a more general, integrated approach 
than that presented in Reference [115], but retaining all the basic 
features of the micromechanics and macromechanics approaches, is needed 
for the prediction of initiation and growth of damage in composites. 

The main objective of the present work was to develop such an 
integrated approach. 
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SECTION 3 


DEVELOPMENT OF AN INTEGRATED FRACTURE CRITERION 


A method of developing a fracture criterion for heterogeneous 
anisotropic composites has been suggested by Wu [124], by re-examining 
the Griffith energy balance equation for homogeneous isotropic mater- 
ials. By expressing the energy terms with greater generality, the 
fundamental assumptions and constraints of classical fracture mechanics 
can be relaxed to allow for material heterogeneity and anisotropy. 

Griffith's energy balance equation for the instability of a crack 
in a brittle material is [1,2] 


dW _ dU 
dA dA - y 


(9) 


where W is the potential energy of the external forces, U is the elastic 
strain energy, y is the surface energy per unit area of the crack 
surfaces and dA is the crack extension. 

Irwin and Orowan [3, A] later modified the Griffith energy balance 
equation to account for the plastic deformation that occurs near the 
crack tip in most engineering materials as follows: 


dW __ dU dU^_ 
dA dA - dA U 


( 10 ) 


where U* is the irreversible strain energy due to plastic deformation. 
Equation (10) is the basis of classical linear elastic fracture mech- 
anics (LEFM) . 



The left-hand side of the Eq. (10) is the input energy rate that is 
released during an incremental crack growth, and the right-hand side is 
the energy absorption rate during the crack extension. The left-hand 
side of Eq. (10) is a function of loading conditions, geometry of the 
body and the orientation of the crack, whereas the right-hand side is a 
constant for a given material. 

The underlying principle in the above mentioned criterion is that a 
crack in a continuum tends to extend when the left-hand side of Eq. (10) 
(i.e., the energy release rate G) reaches a critical right-hand side 
value (G ) , which is a material property, i.e., when 

L 


G = G c (11) 

The above condition, i.e., Eq. (10), which has been derived from 
energy principles, is independent of the constitutive properties of the 
continuum. Therefore, the application of a general fracture criterion to 
composites requires only that the energy terms in Eq. (10) be redefined 
so as to account for the material heterogeneity and anisotropy. 

To begin, we will assume that the composite is a homogeneous, 
anisotropic continuum, containing uniformly distributed microscopic 
flaws, the size and density of which are characteristics of the material 
and manufacturing process. In the vicinity of a microflaw tip (a high 
energy region), we assume that a number of microcracks initiate and 
propagate steadily with increasing external load. The microcracks 
subsequently coalesce and grow in an unstable manner, leading to an 
extension of the macrocrack. 
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The microscopic flaws in a reinforced composite are typically in 
the form of (see Figure 2) : 

a) broken fibers 

b) matrix cracks 

c) fiber-matrix interface debonds 

d) weak sites in the fiber. 

As mentioned earlier, the size and density of the different types of 
flaws depend on the constitutive materials and the manufacturing pro- 
cess. Since the size and density of the microflaws significantly 
influences the failure process in composites, it is necessary to include 
them as material parameters in a fracture criterion for these materials. 

Assuming that the fibers in a unidirectional lamina are packed in a 
regular square array, we can represent all different types of flaws as 
shown in Figure 2. Because of the assumed periodicity, it is possible to 
isolate a repeating cell, such as the one containing a single broken 
fiber surrounded by an array of continuous fibers. A length of four 
times the fiber diameter has been shown [71] to be appropriate for the 
repeating cell in the fiber direction for complete load transfer from 
the broken fiber to the surrounding matrix. 

The sizes of the microflaws within a repeating cell can be non- 
dimensionalized with respect to their individual maximum attainable 
values inside the cell. The microflaw parameters are defined as follows: 
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Figure 2. Typical Microscopic Flaws in a Unidirectional Lamina with Fibers Packed in a 
Square Array. 




a 

mcpr 


Surface area of the matrix crack perpendicular to the fibers 
Maximum attainable size inside the repeating cell 


a 


fmid 


Surface area of the fiber-matrix debond 

Maximum attainable size inside the repeating cell 


The above parameters, along with the fiber volume content (V f ) and 
the mechanical properties of the constituent materials, uniquely define 
the unidirectional composite. 

The repeating cell, which is a representative region of the uni- 
directional lamina, is an ideal model for a micromechanics analysis. But 
the two-material configuration, the presence of microflaws, and the 
complex boundary conditions can pose serious difficulties for a 
continuum analysis. All the above complexities can be easily incor- 
porated in a numerical analysis, however, using a three-dimensional 
finite element approach. In such a numerical analysis it is also pos- 
sible to simulate the onset and growth of microcracks from the micro- 
flaws. The energy absorption during the crack propagation can then be 
evaluated . 

In accordance with our earlier assumption that the energy absorp- 
tion during a macrocrack growth in a composite is in the form of fiber 
breaks, matrix cracks and fiber-matrix interface debonds instead of only 
surface energy of the crack surfaces as in the case of homogeneous 
isotropic materials, the term p in Eq. (10) should be replaced by y, 
which is defined as 

^ dA^frbk + U mopr + U mcpl + U fmid^ 

where = energy absorption due to fiber breaks 

^mcpr = ener Sy absorption due to matrix cracks perpendicular to 
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the fiber direction 


U 


mcpl 


U fmid 


= energy absorption due to matrix cracks parallel to the 
fiber direction 

= energy absorption due to fiber-matrix interface debonding 


The energy balance equation for crack instability in a fiber rein- 
forced composite can thus be written as 


dW 

dA 


dU dlV_ d_ 
dA - dA dA 


< U frbk + 


0 


mcpr 


+ U 


mcpl 


+ U fmid> 


(14) 


The left-hand side of Eq. (14) is the global elastic energy release 
rate, which can be evaluated from a macromechanics analysis, treating 
the composite as a homogeneous anisotropic continuum. The right-hand 
side of the equation is the global energy absorption rate, which is a 
material property. It depends on the size and density of the microflaws, 
and the constituent material properties. The right-hand side should be 
evaluated by means of a rigorous micromechanics analysis, preserving the 
material heterogeneity. Eq. (14) therefore forms the basis of an 
integrated micromechanical and macromechanical fracture criterion 
(IMMFC) for composites. 

The IMMFC can be stated as the onset of a macrocrack in a composite 
laminate that occurs when the energy release rate during a virtual crack 
growth, evaluated from a macromechanics analysis while treating the com- 
posite laminate as a homogeneous continuum, reaches a critical value. We 
define the critical energy release rate as that corresponding to the 
onset of an unstable crack growth from a microflaw. 
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In the micromechanics analysis, in which a repeating cell contain- 
ing microflaws is analyzed for the initiation and propagation of micro- 
cracks under fixed boundary displacement conditions, microcracks are 
likely to grow in a slow and stable manner initially up to a certain 
distance, beyond which they become unstable and grow rapidly across the 
model. The point of crack instability can be established from crack 
growth resistance curves (R-curves) . The energy release rate corres- 
ponding to the point of crack instability is a measure of the critical 
energy release rate (toughness) of the material in that particular 
fracture mode. 

Since the critical energy release rate depends on the direction of 
the loading on the micromechanics model, the three-dimensional micro- 
mechanics analysis will yield six different values of critical energy 
release rate, corresponding to the six independent applied stress 
components, viz., G cll , G^, G^, G^, G cl3 and G 12> where sub- 
scripts 1, 2 and 3 denote the material coordinates (Figure 2). These 
critical energy release rates, which have the dimensions of force, can 
be transformed to any other coordinate system by means of a trans- 
formation. 

The onset of the macrocrack will be in the direction in which the 
energy release rate reaches a critical value in that direction. Thus, 
the IMMFC can also predict the direction of fracture in the composite. 
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SECTION 4 


CRACK GROWTH SIMULATION IN A THREE-DIMENSIONAL 
FINITE ELEMENT ANALYSIS 


The application of the integrated micromechanical and macromechan- 
ical fracture criterion (IMMFC) to composites requires an accurate 
estimation of energy absorption during crack propagation from micro- 
flaws, in both the micromechanics and macromechanics analyses. As 
discussed in the previous section, both analyses are three-dimensional 
in nature, and thus require a three-dimensional analysis. 

The finite element method has been previously used in both two- 
dimensional and three-dimensional fracture mechanics analyses. There are 
a number of special elements to represent the stress singularity that 
exists near a crack tip. A review of the application of the finite 
element method to fracture mechanics is given by Gallaghar [131]. 
Fundamental to the study of fracture mechanics is the study of the onset 
and growth of cracks in a continuum. It is possible in a finite element 
analysis to propagate a crack after initiation. Crack propagation in a 
finite element analysis was introduced by Anderson [132]. In his crack 
propagation simulation, Anderson proposed a method of nodal relaxation 
in which the cohesive crack tip forces at node points ahead of the crack 
are cancelled by applying equal and opposite forces to represent the 
crack extension. Anderson used a constant crack opening angle as a 
fracture criterion and was able to simulate stable crack growth in 
elastic-plastic materials. Light and Luxmoore [133] used the energy 
release rate criterion to predict crack growth. The numerical errors 
associated with a finite element simulation of crack growth in an 



elastic-plastic material have been investigated by Bleackley and Lux- 
moore [134]. A "stiffness derivative" finite element technique has been 
developed by Parks [135] for determining elastic crack tip stress 
intensity factors, which has been further used to derive a virtual crack 
extension analysis for elastic materials by Parks [136] and Hellen 
[137], and for nonlinear materials by Parks [138]. Younan, et al. 
[139,140] used a technique similar to the stiffness derivative method 
with the critical energy release rate criterion to study the crack 
propagation in a heterogeneous anisotropic weldment. 

Existing analytical modeling of crack propagation in composite 
materials is very limited. The crack-closure method [81,83] and the 
failed element [61-64] approaches used for delamination crack growth and 
micromechanics analyses, respectively, were already discussed in Section 
2. In the following sections, a crack propagation simulation in a 
three-dimensional, elastoplastic , finite element analysis will be 
presented . 

4 . 1 General Requirements 

Crack growth simulation in a three-dimensional finite element 
analysis is extremely complex. There are nine possible modes in which a 
stationary crack can advance. Use of failure strength theories (based on 
stress or strain fields) has to be totally ruled out as they cannot 
effectively predict the mode of fracture. One cannot resort to the 
failed element approach [61-64] since this will be grossly inaccurate 
for the present analysis. More general criteria, e.g. , the J-integral 
criterion or the strain energy density factor criterion, which have been 
shown to be successful in planar configurations, have practical limita- 
tions in a three-dimensional application. Since both the micromechanics 
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and macromechanics analyses are envisaged to include elastoplastic 
materials, it is not possible to continuously refine the finite element 
grid ahead of the propagating crack as in the nodal grafting technique 
[141] since in an elastoplastic analysis the stress history at fixed 
points in the material must be maintained. This would restrict the crack 
propagation to occur only between the element surfaces of the original 
grid. To allow sufficient freedom to the propagating crack, the finite 
element grid would have to be very fine and uniform. 

The virtual crack extension method [135,136] combined with the 
energy release rate criterion seems to be the most appropriate approach 
to crack growth simulation in a three-dimensional finite element analy- 
sis. However, for elastoplastic materials the energy release rate during 
an increment of crack growth has to be redefined. Turner [22] has shown 
that the energy release rate in the presence of plasticity is higher 
than that estimated from the linear elastic approximation. 

4 . 2 Energy Release Rate in the Presence of Plasticity 

A simple estimation of the elastic energy release rate in the 
presence of plasticity can be made following the general approach of 
Turner [22], which is a direct extension of the well-known compliance 
method of calculating energy release rate G in the LEFM theory. Refer- 
ring to Figure 3, the path OAB represents a typical nonlinear load- 
deformation curve for an edge cracked elastoplastic plate. If the 
effects of hysteresis are ignored, the elastic unloading path will be 
parallel to the initial loading path OA. If instead, an incremental 
crack growth Aa occurs at constant displacement, in elasticity theory 
the load drops to D due to a reduction in area tAa, where t is the 

thickness of the plate, and BD = -a tAa, where a is the net section 
r n n 
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Figure 3. Estimation of the Elastic Strain Energy Release 
Rate I for an Elastoplastic Crack Growth. 
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stress on the cracked plane, and rises back to D* due to an increase in 
the net stress, i.e., DD' = (w-a)tdo n> where w is the width of the 
plate. The ordinates BD and DD' can also be explained as follows. At 
point B, force equilibrium requires that 

P = a n (w-a)t (15) 

where P is the remote applied force. 

When the crack growth occurs under a condition of zero boundary 
displacement, the crack length, force, and net section stress change to 
a+Aa, P', and o^+do^ respectively. Then force equilibrium requires that 

P* = (a n +da n ) [w-(a+Aa) ] t (16) 

The change in the force is 

P'-P = (a + da ) [w - (a + Aa)]t - a (w - a)t 
n n n 

which can be rewritten as 

P-P' = a tAa - da (w - a)t + da tAa (17) 

n n n 

Neglecting the second order term da R tAa Eq. (17) reduces to 

P-P' = a tAa - da (w - a)t (18) 

n n 
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in which a tAa is the decrease in the force due to the incremental crack 
n 

growth, and da^Cw - a)t is the increase in the force due to the increase 
in the net section stress. 

If the material is assumed to behave elastic-perf ectly plastic, the 

net stress a will be restricted to the yield stress a and there will 
n y 

not be an increase in the force, i.e., DD’ = 0. Unloading from point D 
will be parallel to the increased elastic compliance line D’C. In the 
case of perfect plasticity, the unloading line from D does not pass 
through the transposed origin C. The elastic strain energy release rate 
in the presence of plasticity is represented by area BCED. The area BCD* 
corresponds to the strain energy release rate if the crack growth is 
purely an elastic event. The energy release rate corresponding to the 
area D’CED is the additional energy available in an elastoplastic crack 
growth. It is also possible to actually include the strain hardening of 
the material in the plastic zone near the crack tip instead of assuming 
perfectly plastic response, in which case there will be an increase in 
the net stress, reducing the total energy release rate. For all prac- 
tical purposes, this reduction in the area D’CED due to strain hardening 
will be very small and can be neglected. The elastic strain energy 
release rate in the presence of plasticity as defined here has been 
denoted by I by Turner [22], to distinguish it from the energy release 
rate in the LEFM theory. 

4.3 Fracture Modes and the Estimation of Local Energy Release Rates 

The loading conditions corresponding to the three fundamental 
fracture modes near a crack front are shown in Figure 4. The nine 
possible modes in which a crack in a three-dimensional finite element 
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MODE I 


MODE II 


MODE III 



Figure 4. The Three Fundamental Modes of Fracture Near the 

Crack Tip in a Notched Composite Laminate as Represented 
in a Three-Dimensional Finite Element Model. 
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model can advance are shown in Figures 5 through 7. The nine modes of 
fracture are designated as 


^II 




“yi* “yii* “yiii* m zi’ M ZII and M ZIII* 


As noted in the previous section, we can estimate the energy 
release rate in the presence of plasticity from the rate of change of 
compliance during an incremental crack growth. The relation between the 
rate of change of compliance and the energy release rate during an 
incremental crack growth in a linear elastic material is given by 


G = 


L__ P 2 dC 
2t da 


(19) 


where P is the applied force and C is the compliance of the body (the 
inverse of the stiffness K). 

If is the compliance of the body before an incremental crack 
growth and is the compliance after an incremental crack growth of an 
amount Aa, then the rate of change of compliance can be written as 



da Aa 


Thus, 


or 


C-ip 2 



tAa 




tAa 



( 20 ) 

( 21 ) 

( 22 ) 


where and ^ are stiffnesses corresponding to and C,^. 

Computation of energy release rate in a finite element analysis 
directly using Eq. (22) will involve the solution of the global system 
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Figure 5. Modes M , M v and M with Virtual Crack Plane Normal to the X Coordinate Direction. 
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Figure 6. 


Modes M^, and M Ym wlth Virtual Crack Plane Normal 


to the Y Coordinate Direction. 





Figure 7. Modes M , M an< ^ 

£j 1. ti 11 * 


with Virtual Crack Plane Normal to the Z Coordinate Direction. 







of equations twice (once before an incremental crack growth and once 
after introducing the crack growth). This procedure is not computation- 
ally practical for a continuous crack growth simulation. 

For the simulation of crack growth in a finite element analysis, 
Parks [135] and Hellen [137] have shown that the crack extension affects 
only a few elements near the crack tip. The virtual crack extension 
method given by Hellen is as follows. The total potential energy n of a 
continuum containing an initial crack of length a, in terms of the 
global stiffness matrix [K] , global displacement vector { q} , and global 
load vector {P} in a finite element formulation is given by 

H = | {q} T [K]{q} - {q} T {P} (23) 

If the crack is extended virtually by an amount 6a, the variation of the 
potential energy is 


6JI = |{q} T [6K]{q} + {6q} T [K]{q} - { 5q} T {P) - {q} T {6P} (24) 

where 

[6K] is the change in the global stiffness matrix 

{ 6q } is the change in the global displacement vector 

{6P} is the change in the global force vector 

If the loading forces are far from the crack and are kept constant 

during the virtual crack extension, 6P = 0. Since [K]{q} = {P}, Eq. (24) 
reduces to 

6n = |{q} T [6K]{q} (25) 
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Then the energy release rate for a plate of thickness t is 

G - - £ t - - Hq} ( 26 ) 

This method is sometimes termed the "stiffness derivative method". Only 
a small portion of the stiffness matrix [6K] is populated since the 
crack extension affects only a few elements near the crack tip. 

A similar observation can be made for crack growth under constant 
boundary displacement (fixed grip) conditions. Since the work done by 
the external forces is zero, Eq. (24) becomes 

6 II = SU s = i{q} T [6K]{q} + {<5q} T [K]{q} (27) 

or 

6n = -|{q} T [6K]{q} + {6q} T {P} (28) 

T 

where U is the strain energy of the system. But {6q} {P} = 0 since the 
s 

displacements of the loading boundary are held constant. 

Making use of the above observations, the computation of the energy 
release rate during an incremental crack growth in a finite element 
analysis can be made in terms of the local stiffnesses and forces at the 
crack tip node points. 

In the present analysis, the energy release rate in the presence of 
plasticity I at the crack tip node points is evaluated using the rate of 
change of compliance method. The energy release rate I is defined in 
terms of the local forces and compliance changes. The local forces 
referred to here are the sum of the corresponding incremental force 
components required to move the crack tip node point through the corre- 
sponding incremental displacement of the node point. The method is 
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similar to that used in evaluating the energy release rate by the crack 
closure integral method [85]. In Figure 8, typical force-displacement 
relations corresponding to the three fracture modes (as indicated in 
Figure 4) are shown for an elastoplastic material. 

In the incremental displacement finite element analysis, the linear 
elastic total stiffness K and reduced stiffness K required in calcu- 
lating the local energy release rate can be determined in the first load 
increment, by assembling the contribution of the stiffness coefficients 
of all elements sharing a node point, before and after incremental crack 
growth, respectively. 

The values of I for the three modes of fracture (Figure 8) are 
given by 


I = -i-[F 2 /K - (F - a y AA) 2 /K ] 

I 2AA 1 z ' zz v z zz 


ZZ' 


I = — i— [F 2 /K - (F - T y AA) 2 /K ] 

II 2AA 1 x ' zx v x zx zx 


(29) 


I = — i— [F 2 /K - (F - T y AA) 2 /K ] 

III 2AA 1 y zy v y zy zy J 


where AA is the area of the incremental crack, F^, F y and 7 ^ are force 

components, K , K , K are the total stiffness coefficients, K , K zx 
z z zx zy 

and K are the reduced stiffness coefficients, and o^, T y x , x y are the 
zy * 

yield stress values. The total elastic strain energy release rate is 


I = 


I I + I II + I III 


(30) 


Although the finite element technique of evaluating energy release 
rate was actually developed here for a general three-dimensional 
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a) MODE I 


b) MODE II 


c) MODE III 


Figure 8. Force-Displacement 


Relations for Mode I, Mode II, and Mode III Crack Propagation. 


analysis, the method can be explained more easily using a 

two-dimensional example. Shown in Figure 9a is a crack tip Node Point d 

surrounded by four finite elements (1,2,3 and 4). Under the first load 

increment the Point d gets displaced to d' by a displacement vector 

u,(u,,v,). The deformation of the crack tip elements is shown by the 

dotted lines. With displacements of all node points held constant, the 

force vector F,(F ., F ,) required to move the Point d to d' is 
d v xd yd' 

evaluated by setting the displacement components u^ and v^ equal to zero 
in the matrix force-displacement relation for the crack tip elements 


{F} = [K] {q} 


(31) 


where [K] is the stiffness matrix assembled from the four element 
stiffness matrices, {F} is the nodal force vector for the four elements, 
and { q } contains the nodal displacements. The ratio of the force com- 
ponents (F , and F ..) to the corresponding displacement components gives 
xd yd 

the required total stiffness coefficients at Node Point d before the 
crack extension. 

In Figure 9b, an incremental crack growth has been introduced by 
splitting the Node Point d' into d^ and d£. Under this new configura- 
tion, the Node Points d^ and d^ get displaced to d^ and d^ by displace- 
ment vectors u,. and u jo . The force vectors required to move d. and d„ 
dl dz 1 Z 

to d^ and d^ can be computed, by setting u^ and u^ equal to zero in 
the corresponding force-displacement relations, similar to Eq. (31), in 
which the matrices are obtained by assembling only the top or bottom two 
element matrices. The ratios of the force vectors to the displacement 
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vectors give the required reduced stiffness coefficients at Node Point d 
after an elastic increment of crack growth, 

4.4 Evaluation of Accuracy 

The finite element technique developed to compute the energy 
release rate during an incremental crack growth can be tested for 
accuracy by comparing the computed values with the values obtained by 
standard LEFM methods. The example problem to be considered here is a 
centrally cracked isotropic plate subjected to an inplane tensile force 
normal to the crack plane. The other methods to be used for comparison 
purposes are i) the compliance method, ii) the rate of change of total 
potential energy method, and iii) the stress intensity factor method. 


dC 

In the compliance method, the rate of change of compliance is 


calculated using Eq. (20). The compliances and are obtained from 


°i *£ 


and 


C 2 = r 

2 


(32) 


where and P^ are the forces on the boundary and 6 is the applied 
constant displacement. The energy release rate is obtained using Eq. 
( 21 ). 

In the rate of change of total potential energy method, the energy 
release rate is given by 


G = 


Hi ~ Il2 

tAa 


(33) 


where the total potential energies before and after an incremental crack 
growth are given by 
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and 


(34) 


n 


1 





- ft 


2 


where and represent the strain energy of the plate before and 
after the crack growth, respectively, and ft^ and are th e correspond- 
ing potential energies of the external forces. 

The finite element grid used to model the centrally cracked plate 
is shown in Figure 10, The plate was assumed to be made of a linearly 
elastic material having a Young’s modulus E of 69 GPa (10.2 Msi) and a 
Poisson’s ratio v of 0,3. The plate dimensions were 15.2 cm (6 in) in 
length, 10.2 cm (4 in) in width, and 0.13 cm (0.05 in) in thickness. The 
two ratios of crack length a to width b considered were 0.2 and 0.3. A 
total of twenty eight 20-node quadratic isoparametric brick elements 
were used. In order to represent the l//r stress singularity near the 
crack tip, the 20-node elements were collapsed into 15-node quarter 
point, triangular wedge singular elements [142]. Tensile stresses were 
applied to the plate by prescribing uniform displacements at the boun- 
dary nodes. As can be seen from the results presented in Table 1, the 
agreement of the present analysis with existing analyses is excellent. 

4 . 5 Crack Growth Simulation 

A three-dimensional, elastoplastic , generally orthotropic finite 
element computer program patterned after that of Monib and Adams [144] 
was developed for use in the present work. The present program includes 
20-node quadratic isoparametric brick elements and a crack propagation 
capability. The basic elastoplastic analysis theory and the structure of 
the present computer code are given in Appendix A. 

At the node points near an existing crack, the energy release rates 
for all possible modes of fracture are evaluated at each load increment, 
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SINGULAR ELEMENT 


Figure 10. Finite Element Grid Used to Model the Centrally Cracked 
Isotropic Plate. 



Table 1 


Comparison of Energy Release Rates During an Incremental 
Crack Growth Obtained by Different Analytical Methods for a 
Central Cracked Isotropic Plate 


Method 


Energy Release Rate 
J/m (lb/ in) 


Compliance Method (Eq. 21) 


1 P 

G = =r ~{C~C 0 ) 

2 Aat 1 2 


6089 34.4 


Rate of Change of Potential Energy (Eq. 33) 


G = 



tAa 


5947 


33.6 


Stress Intensity Factor [143] 

wK t yp/T 

G = — — , where K = — — 6106 34.5 

E I wt 

and Y = 1.82 

Present Analysis 6053 34.2 
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and compared with the corresponding critical energy release rate values 
for the material. When the computed energy release rate at a node point 
is equal to the critical value in an increment, that node point is 
separated into two node points and the reaction forces are applied to 
the old and new node points. The element nodal connectivity of all 
elements sharing the new node point is changed to simulate the crack 
growth in the model. If the crack tip node point at which the energy 
release rate reaches a critical value happens to be a corner node, the 
crack is assumed to extend up to the nearest mid-side nodes; if the node 
point is a mid-side node, the extended crack surface is assumed to be 
that connecting the adjacent corner nodes and the mid-side node on the 
opposite side. The possible forms of crack extension are shown in Figure 
11. Before applying the next load increment, the system of equations 
with new boundary conditions and reaction forces is solved again to 
check for any further crack growth. The crack growth simulation 
technique and its implementation in a three-dimensional elastoplastic 
finite element computer code is further explained in Appendix B. 

4 . 6 Example Problem 

The accuracy of the crack growth simulation technique developed in 
the present analysis was verified by applying it to predict self-similar 
crack growth in a centrally cracked 2024-T3 aluminum plate subjected to 
inplane Mode I loading. 

The finite element grid and plate geometry were assumed as shown in 
Figure 10. The material properties used for the 2024-T3 aluminum alloy 
[143] were as follows: Young’s modulus = 70.3 GPa (10.2 Msi) , Poisson’s 
ratio = 0.345, yield strength = 345 MPa (50 ksi) , and tensile strength = 
485 MPa (70 ksi). The nonlinear properties of 2024-T3 aluminum [143] 
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were incorporated in the analysis by expressing the stress-strain 
relation of the material in terms of Richard-Blacklock curve-fitting 
parameters (see Appendix). The displacement constraints on the collapsed 
nodes of the singular elements were removed, to represent the 1/r type 
of singularity of a perfectly plastic material. Such perfectly plastic 
singular elements are extensively used in elastoplastic fracture 
mechanics [145]. The thickness of the plate was taken as small as 
possible, viz., 0.0025 cm (0.001 in), to simulate plane stress condi- 
tions without introducing errors due to resulting large aspect ratios of 
the grid elements. 

In Figure 12, the predicted crack extension at various applied 
stress levels and the associated extent of plastic yielding are shown. 
The shapes of the plastic yield zones were obtained by drawing smooth 
curves around the Gaussian integration points at which the material was 
predicted to have yielded. It should be mentioned here that in earlier 
trial runs in which greater plate thickness values were considered, the 
predicted yield zone sizes were smaller than those shown in Figure 12, 
and they varied in size through the thickness of the plate. The extent 
of plastic yielding at various applied stress levels predicted in the 
present analysis is in general agreement with those given in Reference 
[146] for a central cracked aluminum plate. 

During the process of crack propagation (particularly for large 
scale yielding) , the analysis also predicted that some of the elements 
which had previously plastically deformed start unloading as the crack 
passes by, as one might expect. 

It was observed that the crack propagated in small increments 
initially, up to about half the plate width, and then propagated rapidly 
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h — n 10,2 mm (0,4 in) K s — H 12.2 mm (0,48 in) ^ ^20,3 mm (0,80 in) 

a) a = 58 MPa (8.5 ksi) b) a = 81 MPa (11.9 ksi) c) a = 199 MPa (28.6 ksi) 

Figure 12. Crack Propagation and Associated Plastic Yielding in a Centrally Cracked 
2024-T3 Aluminum Plate as Predicted by the Present Analysis. 






across the plate. This may have been due to the coarseness of the finite 
element grid used, especially towards the later stages of crack growth, 
but the applied stress a = 330 MPa (47.9 ksi) at the point of observed 
instability of the crack agreed with the critical stress value a ■ = 348 
MPa (50.4 ksi) for the 2024-T3 aluminum [143]. The fracture toughness 
value of 41.8 MPa/m" (38 ksi/In) computed using = /GE, where G is the 
energy release rate at the onset of unstable crack growth and E is the 
Young's modulus of the material, was also in good agreement with the 
material toughness value of 44 MPaVm" (40 ksi/in) reported in Reference 
[143]. 

The predicted crack opening displacements (crack shapes) of the 
propagating crack are plotted in Figure 13. The effect of plastic 
yielding at higher applied stress levels in blunting the crack tip can 
be seen. Neither experimental nor analytical crack opening displacement 
values were available in the literature for the material and crack 
geometry used here, preventing comparisons with the crack opening dis- 
placements predicted here. 


56 
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Figure 13. Crack Opening Displacements in the Centrally Cracked 

2024-T3 Aluminum Plate at Various Applied Stress Levels. 
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SECTION 5 


MICROMECHANICS ANALYSIS 


5.1 Model Geometry 

Assuming that the fibers in a unidirectional lamina are packed in a 
regular square array, a repeating cell can be identified. One such is 
that containing a single broken fiber surrounded by continuous unbroken 
fibers, as previously indicated in Figure 2. In this configuration, the 
density of the broken fibers is 25 percent of the total number of 
fibers. Other types of microflaws to be considered include matrix cracks 
parallel to the fibers, and fiber-matrix interface debonds. These 

microcracks are assumed to be distributed symmetrically with respect to 

three orthogonal planes. 

Figures 14 and 15 show the two three-dimensional finite element 
models developed for the analysis. Taking advantage of the assumed 

symmetry, only one octant of the repeating cell need be modeled. The 

first model (Figure 14) contains a total of 282 20-node, quadratic 
isoparametric elements, whereas the second model (Figure 15) contains 
144 similar elements. Due to the memory (including extended core) limits 
of the CDC Cyber 760 computer available for use at the University of 
Wyoming, it was not possible to run the first, more refined model. An 
attempt was made to run this model on the Department of Mechanical 
Engineering’s Prime 550 minicomputer, but the running time per increment 
was prohibitively long. However, the limited results for the few incre- 
mental solutions which were obtained using the Prime 550 (with double 
precision) were later used to evaluate the much coarser second model 
(Figure 15). Stresses predicted using the two different models differed 



1 (Z) 



2 (X) 


Figure 14. Finite Element Grid Developed for the Micromechanics 
Analysis. 
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Figure 15. Finite Element Grid with Representative Microflaws Used 
in the Present Micromechanics Analysis. 
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by only 4 percent (maximum). The very large aspect ratios of the inter- 
face elements and the distortions of the fiber elements in the second 
model did not introduce significant errors. 

There are 945 node points in the second model (1583 in the first). 
Three different microflaws, in the form of broken fibers, matrix cracks 
parallel to the fibers, and fiber-matrix interface debonds, have been 
represented in this model, by altering the boundary conditions. The 
material microflaw parameters (defined by Eqs. 12 in Section 3) used in 
the present analysis were: 


a f rbk 
a mcpl 
a fmid 


0.250 

0.020 

0.017 


The model incorporates a double-node concept at the junction of the 
broken fiber and the surrounding matrix, in order to represent the 
actual conditions of discontinuity of the fiber at the break, while 
retaining the continuity of the matrix material at the same point. The 
double nodes are constrained to have identical displacements in the x 
and in the y directions. 

The displacement boundary conditions for the octant of the repeat- 
ing cell shown in Figure 15 are that the normal displacements of all 
nodes on the vertical boundaries be uniform (from symmetry), and that 
the normal displacements of all nodes on the upper horizontal boundaries 
(surface EFGH) also be uniform (generalized plane strain loading). This 
complex boundary condition has been achieved in the present three- 
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dimensional model by assigning the same global degree of freedom to all 
nodes that are required to have same boundary displacements. 

5. 2 Material Properties 

The constituent materials used for the unidirectional lamina were 
Hercules AS4 graphite fibers and Hercules 3501-6 epoxy matrix. The 
temperature- and moisture-dependent properties of the neat (unrein- 
forced) matrix material as obtained by means of solid rod torsion tests 
[147] have been expressed in the present analysis in a polynomial form 
as 


P(T,M) = CjT + C 2 M + C 3 TM + C 4 (35) 

where P is the material property of interest, T is temperature (°C), M 
is the moisture content (weight percent) and the C's are coefficients of 
the polynomial (see Table 2). Temperature- and moisture-dependent 
octahedral shear stress-octahedral shear strain curves for the Hercules 
3501-6 resin are shown in Figure 16. The matrix is treated as an iso- 
tropic, elastoplastic material in the analysis, by using a three- 
parameter Richard-Blacklock model [148]. 

The properties of the graphite fibers used in the analysis are 
given in Table 3. The anisotropic nature of the fibers has been taken 
into account in the present analysis, as explained in the Appendix. 

5 . 3 Loading Conditions 

The main objective of the present micromechanics analysis was to 

evaluate the six critical energy release rate values, viz, G . G 

’ Cll C22 

^C33 * ^c23’ ^C31 an< ^ ^C12* These critical energy release rate values 
correspond to the loading conditions denoted by their subscripts. In 
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Table 2 


Coefficients of the Polynomial (Eq. 35) Used to Define the 
Temperature-and Moisture-Dependent Properties of the Hercules 

3501-6 Epoxy Resin. 


Property = CjT + C 2 M + C 3 TM + C 4 


Property C^ 


Shear Modulus, G —1.14 x 10 3 -2.49 x 10 3 -3.38 x 10^ 2.60 x 10 3 

(psi)* 

Richard 

Blacklock 

Curve-Fit „ 

Parameters n* 2.56 x 10 2.78 x 10 1.65 x 10 1.54 x 10 

to (psi)* -1.50 x 10 2 -1.36 x 10 3 6.46 x 10° 2.67 x 10 4 

Ultimate Shear -8.73 x 10^ -6.03 x 10 2 1.86 x 10^ 1.66 x 10 4 

Strength, x^^Cpsi)* 

Coefficient of 1.22 x 10~ 7 1.04 x 10 -6 -5.90 x 10 -10 3.83 x 10" 

Thermal Expansion, 
a ( # C) 

Coefficient of 0 0 0 3.20 x 10 

Moisture Expansion, 

8 (%M) _i 

Poisson's Ratio, v 0 0 0 0.34 x 10*^ 


*Based on solid rod torsion shear data 
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Figure 16 . Temperature— and Moisture— Dependent Octahedral Shear Stress— Octahedral Shear 
Strain Curves for Hercules 3501-6 Epoxy Resin, as Obtained from Solid Rod 
Torsion Tests [ 147 ]. 
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Table 3 


Mechanical Properties of Hercules AS4 Graphite Fiber [147] 


Longitudinal 


Modulus, E^ 

221 GPa 

(32.0 Msi) 

Transverse 
Modulus, E t 

13.8 GPa 

(2.0 Msi) 

Longitudinal Shear 
Modulus, 

34.5 GPa 

(5.0 Msi) 

Transverse Shear 
Modulus, G 

5.52 GPa 

(0.80 Msi) 

Major Poisson’s 
Ratio, v lt 

0 . 

20 

In-Plane Poisson’s 
Ratio, v 

0 . 

25 

Longitudinal Tensile 
Strength, 

3100 GPa 

(450 ksi) 

Transverse Tensile 
Strength, a 

345 GPa 

(50 ksi) 

Longitudinal Shear 
Strength, x£ 

1550 GPa 

(225 ksi) 

Transverse Shear 
Strength, x 

172 GPa 

(25 ksi) 


Longitudinal Coefficient 

of Thermal Expansion, -0.36x10 /°C 


Transverse Coefficient 
of Thermal Expansion, 


18. Ox 10~ 6 /°C 
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general, to evaluate all six values, the micromechanics model must be 
subjected to six separate loading conditions, viz, 0 ^, °22’ °33’ T 23’ 
and T i2* however, since the model and the microflaws are assumed to 
be symmetric about the x-z and y-z planes (in the transverse direction), 

G C22 = G C33 and G C13 = G C21* 

This reduces the total number of loading conditions to four (viz. 
Off, a 22 , t and t^)- Usin S the finite element model of the octant of 
the repeating cell (Figure 15) , the longitudinal and transverse normal 
loads (see Figures 17a and 17b, respectively) are easily applied as 
prescribed boundary displacements, in increments, whereas the two 
longitudinal shear loading conditions (see Figures 17c and 17d) pose 
difficulties. Under shear loading, the assumption of three orthogonal 
planes of symmetry used in modeling an octant of the repeating cell is 
not valid. Longitudinal shear loading (Figure 17c) can be modeled dir- 
ectly, however, as presented in detail in Reference [40]. The longitud- 
inal shear loading is applied by prescribing uniform displacements to 
all node points on face FBCG (Figure 15 or Figure 17c) in the longi- 
tudinal direction and removing the uniform displacement constraint on 
all node points on the top face (EFGH) except along lines EH and FG. A 
general approach to apply transverse shear loading in a two-dimensional 
micromechanics analysis is to identify a repeating cell with its sides 
at 45° to the transverse co-ordinates (2 and 3 in Figure 15), and then 
to apply equal and opposite normal stresses to the cell boundaries. It 
is not possible to identify such a repeating cell in the present case 
that has same percentage of broken fibers. However, the crack propaga- 
tion due to shear loading in the transverse plane has been studied in 
the present analysis by applying equal and opposite normal stresses in 
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the two transverse directions (Directions 2 and 3 in Figure 17d) . This 
is equivalent to transverse shear loading at 45° to the transverse 
co-ordinates • 

Before applying any mechanical loads, the thermal conditions 
induced after curing the graphite/epoxy composite at (177°C) and then 
cooling down to room temperature (21°C) were simulated by reducing the 
stress-free reference temperature of 177°C to room temperature in 10 
equal increments. 

5 . 4 Crack Initiation and Propagation 

For the three different types of microflaws considered in the 
present analysis (broken fibers, matrix cracks parallel to fibers, and 
fiber-matrix debonds) it could be judged beforehand that each individual 
flaw type would become critical only under certain loading conditions 
(e.g., the broken fiber will become critical only under longitudinal 
tension and longitudinal shear loading conditions) • These observations 
considerably reduced the total number of computer runs required. The 
different combinations of types of microflaws and loading conditions 
actually considered are summarized in Table 4. 

Table 4 

Loading Conditions Considered with Different Types 
of Microflaws in the Present Analysis 

Flaw Type Loading Conditions 

Broken Fibers o ^ “ 

Matrix Cracks - o ^ - t ^ 

Fiber-Matrix Debond - o ^ - T 23 
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5.4.1 Broken Fibers 


The fiber break is introduced in the finite element model (Figure 
15) by simply removing the displacement constraint of all node points in 
the fiber lying in the plane of the break. If the fiber break is intro- 
duced before the cooling temperature increments are applied, the graph- 
ite fiber, which has a negative longitudinal coefficient of expansion, 
tends to extend beyond the symmetry plane. This has been prevented in 
the present analysis by not introducing the flaw until the system is 
cooled down to room temperature. The crack propagation flag is also 
turned off during the incremental loading. 

The applied longitudinal displacement increments used were 0.0076 
cm (0.003 in). In Figure 18, the predicted crack propagation at various 
applied stress levels is schematically shown. At the applied stress of 
284 MPa (50 ksi) , all nodes near the initial crack opened in modes M z , 
M and M TTT (predominantly in M ), forming a radial crack front. 
During subsequent load increments, the crack steadily grew until it 
reached the diagonally opposite fiber. The crack then spread all across 
the matrix. Although not shown here, plastic yielding was restricted to 
small volumes near the propagating crack. Comparisons of the sizes of 
the plastic yielding zones indicated that as the crack front advanced, 
the nodes previously yielded started unloading. 

In Figure 19, the predicted crack growth from the broken fiber due 
to longitudinal shear loading is shown. The applied longitudinal dis- 
placement increments on the vertical face were 0.0076 cm (0.003 in). Two 
radial cracks initiated in the plane of the fiber break, and propagated 

V 

around the adjacent fibers. 
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a) prior to loading b) 8.3 MPa (1.2 ksi) 



c) 53 MPa (7,8 ksi) d) 61 MPa (9.0 ksi) 

Figure 19. Initiation and Growth of Cracks from a Broken Fiber 
as a Function of Longitudinal Shear Applied Stress. 




5.4.2 Matrix Cracks Parallel to Fibers 


Two matrix cracks parallel to the fiber direction (Figure 20a) were 
introduced into the finite element model prior to loading by splitting 
the central node points on each of the sides into two nodes. The matrix 
cracks are symmetrically oriented with respect to a diagonal vertical 
plane, thus maintaining the transverse symmetry of the model. 

Considering transverse normal loading first (see Table 4) , the 
transverse normal incremental displacements applied were 0.0051 cm 
(0.002 in), on the face normal to the 2-axis (see Figure 20). Figure 20 
shows the predicted crack propagation at various applied stress levels. 
The matrix crack normal to the load direction propagated in its initial 
plane up to the center of the model, and then suddenly changed its 
direction and propagated at 45° to the loading direction. In subsequent 
load increments, the crack growth was in the fiber direction. In an even 
later stage, the other initial matrix crack (parallel to the loading 
direction) started to propagate in its initial plane, all the way across 
the model (see Figures 20c and 20d) , and then in the fiber direction. 

The initiation and growth of microcracks from the initial matrix 
cracks under transverse shear loading is shown in Figure 21. The incre- 
mental displacements on the transverse faces were +0.0051 cm (+0.002 
in) and -0.0051 cm (-0.002 in) in the 2- and 3-directions, respectively. 
As in the case of transverse tension, the initial crack, after running 
in its plane for some distance, changed its direction by 45° to the 
transverse plane. 

The mixed mode crack growth predicted when the unidirectional 
composite containing matrix cracks is subjected to either transverse 
tension or transverse shear is important from another viewpoint also. 
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The crack growth predicted may very well explain the experimentally 
observed zig-zag path of the delamination. It can be concluded that 
cohesive matrix cracks in the "resin-rich** zone at the ply interface are 
influenced by the presence of adjacent fibers, and change their path 
during propagation As a result, the "fracture toughness** of the com- 
posite also increases. This is in agreement with the experimental 
observation that the fracture toughness of the reinforced resin is 
higher than that of the neat resin. The transverse tension analysis can 
be directly correlated with the double-cantilever beam test in which a 
delamination is introduced between plies. 

5.4.3 Fiber-Matrix Debond 

In this flaw model, one of the fibers in the finite element model 
of the repeating cell was partially debonded from the matrix material by 
introducing double node points at the fiber-matrix interface in the x-y 
symmetry plane. The model was then subjected to transverse tension and 
transverse shear loads as in the case of the matrix crack flaw model. 
The crack growth at various stages of transverse tensile loading is 
shown in Figure 22. The crack growth pattern for transverse shear 
loading was similar. Only a part of the initial debond extended in the 
fiber direction, as shown in Figure 22. 

5.5 Critical Energy Release Rates 

The microcracks from the three different types of initial micro- 
flaws considered grew in a number of small increments. By further 
refining the finite element grid and reducing the load increment sizes, 
it should be possible to obtain slow and stable crack growth from 
microflaws. The fact that additional load increments are required to 
propagate a crack suggests that there is increasing resistance to crack 
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a) prior to loading 


b) M MPa (2,1 ksi) 



c) 26 MPa (3.8 ksi) d) 62 MPa (9.1 ksi) 


22. Initiation and Growth of a Crack from an Initial Fibsr 
Matrix Debond as a Function of Transverse Tensile 
Applied Stress. (Fiber on Vertical Axis Removed for 
Viewing Clarity) 



growth. This situation is similar to that observed in the case of thin 
metal plates. The slow and stable crack growth in metals has been 
explained in terms of R-curves [13,149]. The energy release rate 
corresponding to the point of crack instability can be interpreted as 
the critical energy release rate or fracture toughness. However, in the 
present analysis, the loading increments were not small enough to 
generate adequate R-curves. Therefore, the energy release rate 
corresponding to only the first or second increment of crack growth was 
taken as the critical energy release rate here. The energy release rate 
was evaluated by dividing the difference between the global strain 
energy before and after an incremental crack growth by the total area of 
the incremental crack. The potential energy of the external forces was 
zero since the crack growth was modeled as occurring under fixed grip 
conditions. The computed critical energy release rates for the different 
types of microflaws and loading condition combinations are tabulated in 
Table 5. 


Table 5 

Computed Critical Energy Release Rates from the Micromechanics Analysis 


Micro flaw 
Type 

Critical Energy Release Rate, j/m (in-lb/in ) 

G C11 G C22 G C13 G C23 

Broken 

j 

Fibers 

434 (2.45) — 1873 (10.58) 

Matrix 


Cracks 

262 (1.48) -- 273 (1.54) 

Fiber- 


Matrix 


Debond 

- 

119 (0.67) -- 350 (1.98) 
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For use in the macromechanics analysis to be presented in the next 
section, the lower of the values for a particular load increment was 
taken as the critical value. The six critical energy release rates 
actually used in the subsequent macromechanics analysis are given in 
Table 6. 


Table 6 

Computed Values of Critical Energy Release Rates Selected for 
Use in the Macromechanics Analysis 

2 2 
Parameter (J/m ) (in-lb/in ) 


G cn 

434 

2.45 

G C22 

119 

0.67 

G C33 

119 

0.67 

G C23 

273 

1.54 

G C13 

1873 

10.58 

G C12 

1873 

10.58 
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SECTION 6 


MACROMECHANICS ANALYSIS 

6.1 Notched Composite Laminates 

The numerical results of a macromechanics analysis of the initia- 
tion and growth of a delamination crack near a through-the-thickness 
notch crack in a composite laminate is presented in this section. In the 
present macromechanics analysis, the individual laminae are treated as 
homogeneous anisotropic continua. A central notched and a single-edge 
notched [±45/0] graphite/epoxy laminate were considered. It was assumed 
that delamination involved only matrix-dominated fracture, and that the 
crack surface is parallel to the ply interface. 

The matrix-dominated transverse normal, transverse shear and 
longitudinal shear properties of the lamina used in the present analysis 
are shown in Figure 23. These properties were generated using a two- 
dimensional finite element micromechanics analysis [147]. The nonlinear 
stress-strain response of the ply material indicates the need to use the 
elastoplastic stress analysis and fracture criterion for crack growth 
developed in the present analysis. 

6.2 Finite Element Model 

The finite element grid used to model the central and single-edge 
notched [±45/0] s graphite/epoxy laminate is shown in Figure 24. The 
laminates were assumed to be 15.2 cm (6 in) in length, 10.2 cm (4 in) in 
width, and to have a notch length a to plate width b ratio of 0.2 in 
both cases. Twenty-node, quadratic, isoparametric elements were used. 
Each individual ply was modeled by a separate layer of elements. The 
model shown previously in Figure 10 is a single-layer model. Even though 
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a) transverse tension 


b) transverse shear 



c) longitudinal shear 


Figure 23. Matrix-Dominated Properties of a Unidirectional Hercules AS4/3501-6 Graphite/Epoxy 
Composite as Predicted by a Finite Element Micromechanics Analysis [147] (fiber 
volume = 60 percent) • 







SINGULAR ELEMENT 


Figure 24. Finite Element Grid Used to Model the Central and Single 
Edge Notched [±45/0] Graphite/Epoxy Laminate, a/b =0.2 




the assumed material symmetry relative to the x-z plane does not actual- 
ly exist in the case of a center notched laminate, it was assumed here 
that this does not affect the delamination crack growth in the one-half 
laminate modeled. In order to reduce the total number of 20-node ele- 
ments required to accurately represent the stress state near the crack 
tip, the 20-node isoparametric elements were collapsed into 15-node, 
quarter-point, triangular, wedge-like, singular elements in this region 
[ 14 2 ] . These elements, if the collapsed nodes at the crack tip are not 
constrained to have the same displacements, represent a perfectly 
plastic (l/r type) singularity [150]. 

6.3 Crack Initiation and Propagation 

The finite element model of the central notched and single-edge 
notched [±45/0] graphite/epoxy laminate was subjected to inplane 
uniaxial tensile stresses, applied in small increments by prescribing 
boundary displacements. At each increment of loading, the local elastic 
energy release rates were calculated at all node points lying near the 
crack tip in the +45/-45 and -45/0 interfaces. If the computed energy 
release rate at any of these nodes was equal to the critical energy 
release rate, the node was split into two nodes and the reaction forces 
were applied to the split nodes. Holding the boundary node displacements 
constant, the system of equations was solved again to check for further 
crack extension due to the reaction forces at the split nodes. Experi- 
mentally measured [100] fracture toughness values for a Hercules 

2 

AS4/3501-6 graphite/epoxy composite, viz, G = 130 J/m and G = 

IG IIC 

2 

= 230 J/m , were used as critical energy release rates in the 
present analysis. 


84 



Figures 25 and 26 show the growth of the plastic zone near the 
crack tip between +45/-45 and -45/0 interfaces, for the central notched 
laminate during different stages of loading. Figures 27 and 28 show the 
growth of plastic yielding for the single-edge notched laminate. In both 
the central notched and single-edge notched laminates, the plastic 
yielding was confined to a very small region near the through-the- 
thickness crack tip. 

Figures 25 through 28 presented the extent of interface yielding. 
The onset and growth of actual delamination cracks between the +45/-45 
and -45/0 interfaces for the case of the central notched laminate are 
shown in Figures 29 and 30. The applied stress increments were kept 
small enough to cause only one increment of crack growth in an increment 
of applied stress. A particular mode (Mode I, II or III) was predominant 
in each incremental crack advance; this is indicated in Figures 29 and 
30 by different shadings to represent the crack growth. The delamination 
crack growths in the case of a single-edge notched laminate are shown in 
Figures 31 and 32. In both the central and single-edge notched lami- 
nates, the delamination was more extensive between -45/0 interfaces than 
between +45/-45 interfaces. 

6.4 Crack Growth Resistance and Fracture Toughness 

In all cases, for loadings beyond the maximum applied stress levels 
shown in Figures 29 through 32, the analysis predicted a very large area 
of delamination. This was undoubtedly due to the very coarse finite 
element grid used in the present analysis. With a finer grid, the crack 
growth would likely continue further in the same slow and stable manner 
as observed during the initial phase in the present analysis. One of the 
explanations given for such stable crack growth behavior is that the 
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a) 120 MPa (17.4 ksi) 


b) 223 MPa (32.3 ksi) 


c) 291 MPa (42.2 ksi) 


Figure 25. Material Yielding at Increasing Levels of Applied Axial Tensile Stress Near the 
Crack Tip at the +45 / —45 Interface in a Central Notched [ ±45 / 0 ] s Graphite/Epoxy 
Laminate. 




a) 120 MPa (17,4 ksi) b) 223 MPa (32,3 ksi) c) 291 MPa (42.2 ksi) 


Figure 26. Material Yielding at Increasing Levels of Applied Axial Tensile Stress Near the 
Crack Tip at the -45/0 Interface in a Central Notched [±45/01 g Graphite/Epoxy 
Laminate. 




a) 120 MPa (17.0 ksi) 


b) 223 MPa (32.3 ksi) 


c) 291 MPa (02,2 ksi) 


Figure 27. Material Yielding at Increasing Levels of Applied Axial Tensile Stress Hear the 

Crack Tip at the +45/- 45 Interface in a Single-Edge Notched [ ±45/0] g Graphite/Epoxy 
Laminate. 







a) 120 MPa (17,4 ksi) b) 223 MPa (32,3 ksi) c) 291 MPa (42,2 ksi) 


Figure 28. 


Material Yielding at Increasing Levels of Applied Axial Teiisile ohite/Epoxy 

Crack Tip at the -45/0 Interface in a Single-Edge Notched [±45/0] g Graphite/Epoxy 

Laminate. 



a) 120 MPa (17.7 ksi) 


b) 171 MPa (24. B ksi) 
M MODE II 


c) 223 MPa (32.3 ksi) 
M MODE III 


figure 29 . Delamination Growth Near the Crack Tip at Che + *^'^ of t AppUed 1 Axial C Tensile°S tress 
[ ± 45/0] Craphite/Epoxy laminate at Increasing Levels of Applied Axial 





o 



0 



b) 120 MPa (17.4 ksi) 


0 





\ 

// 


/7 

V\ 






c) 171 MPa (24,8 ksi) 


MODE I 


MODE II 


MODE III 


Figure 30. 


Delamination Growth Near the Crack Tip at the -45/0 Interface in a Central Notched 
[±45/0] Graphite/Epoxy Laminate at Increasing Levels of Applied Axial Tensile Stress 
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c) 223 MPa (32,3 ksi) 
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Figure 31. Delamination Growth Near the Crack Tip at the -45/0 Interface in a Central Notched 

[±45/0] Graphite/Epoxy Laminate at Increasing Levels of Applied Axial Tensile Stress, 
s 
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B) 171 MPa (24.8 ksi) 
M MODE II 


a) 120 MPa (17.4 ksi) 
■ MODE I 


c) 223 MPa (32.3 ksi) 


HI MODE III 


Figure 32. Delamination Growth Near the Crack Tip at the +45/-45 Interface in a Single-Edge 
Notched [±45/0] Graphite/Epoxy Laminate at Increasing Levels of Applied Axial 
Tensile Stress. 




energy absorption rate continuously increases with increasing plastic 
zone size in front of the crack tip [95]. In the present analysis, the 
elastic strain energy release rate at the next node point attains the 
critical value only after sufficient material yielding. The delamination 
crack growth observed in the present analysis can be characterized by 
R-curves [150,151], and a point of instability can thus be established. 
The energy release rate corresponding to the onset of unstable crack 
growth is then defined as the delamination fracture toughness of the 
particular notched composite. Owing to the coarseness of the present 
finite element grid, crack growth beyond the singular crack tip elements 
is not considered accurate. Nevertheless, for general information, the 
energy release rates corresponding to the onset of the large delamina- 
tions observed in the present model are given in Table 7. It should be 
noted that in the present analysis, only delamination crack growth was 
allowed; all other possible modes of crack growth, e.g., transverse 
intralaminar cracking, were completely suppressed. 


Table 7 

Critical Energy Release Rates Predicted During Delamination Crack 
Growth in Notched [±45/0] Graphite/Epoxy Composite Laminates 

O 


Notch Geometry 


Critical Energy Release Rate 
(J/m 2 ) (lb/ in) 


Central Notch 
Single-Edge Notch 


515 

2.92 

511 

2.89 
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SECTION 7 


INTEGRATED ANALYSIS 

The ability of the integrated micromechanical and macromechanical 
fracture criterion (IMMFC) and the crack propagation scheme developed in 
the present analysis to predict the onset and growth of cracks in 
general multilayered composite laminates will be demonstrated in this 
section by means of an example problem. A single-edge notched [±45/0] g 
graphite/epoxy laminate is assumed to be subjected to inplane tensile 
stress normal to the notch (as shown previously in Figure 24) . Individ- 
ual laminae of the laminate are assumed to contain uniformly distributed 
microflaws of the same size and density as used in the micromechanics 
analysis of Section 4. In accordance with the IMMFC, the critical energy 
release rate values evaluated by a rigorous micromechanics analysis of 
the unidirectional lamina, as listed previously in Table 6, were used as 
the criterion for the onset of cracking. These G c values, which were 
derived in material coordinates, were transformed here into laminate 
coordinates for the +45 and -45 plies. 

The finite element grid, the layup of the laminate, and the notch 
dimensions were assumed to be the same as these used for the macro- 
mechanics analysis (see Figure 24 of Section 6). 

7.1 Crack Propagation 

The incremental inplane tensile stress was applied by prescribing 
uniform displacements to all appropriate boundary node points. Displace- 
ments of 0.00127 cm (0.0005 in) and 0.00254 cm (0.001 in) were applied 
in the first two increments, respectively, followed by increments of 
0.00762 cm (0.003 in) in subsequent increments. 



In Figures 33 through 35, perspective views of the total crack/ 
damage growth are shown, at three successively increasing levels of 
loading. It must be noted that the sizes of the incremental cracks are a 
direct function of the finite element sizes. In Figure 36, the average 
laminate stress is plotted against the applied normal tensile strain. 
The several horizontal steps in this stress-strain curve correspond to 
the large-scale crack propagation which occurred at higher load levels. 
Beyond applied strains of 6.75 percent there was no increase in the 
stress level. 

Details of the predicted crack/damage growth between and within 
individual plies at various applied stress levels are schematically 
shown in Figures 37 through 39. Due to the complex nature of the pre- 
dicted crack growth, results are shown individually between each ply 
interface, and also on the top surface. It can be seen that, in addition 
to the self-similar growth of the original crack, a number of through- 
the-thickness transverse cracks running parallel to the fibers are 
initiated during the initial stages of loading. The delamination mode 
becomes predominant at the later stages of loading. In the present 
analysis, only node points in a narrow band on either side of the 
original crack were considered for possible cracking. 

No experimental results were available in the literature for the 
specific layup sequence and laminate geometry used in the present 
example problem. However, the predicted strength and crack propagation 
extent are in very good agreement with general experimental observations 
[151]. This example problem demonstrates that the IMMFC and crack 
propagation scheme developed in the present analysis are very effective 
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<W=A> 


SINGLE-EDGE NOTCHED 

C-M5/-45/0I3S GR/EF LAMINATE 
INPLANE TENSION 



Figure 36 . Predicted Stress-Strain Response of a Single-Edge 
Notched [± 45 / 0 ] Graphite/Epoxy Laminate Subjected 
to an Axial Tenlile Stress. 
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STRESS CKSD 





a) between 0/0 PLIES b) between 0/-45 PLIES 


Figure 37. Crack/Delamination Growth in a Single-Edge Notched 

[±45/0] Graphite/Epoxy Laminate at an Applied Axial 
Tensile S Stress of 105 MPa (15.2 ksi) . 
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TRANSVERSE 
CRACKS 
PARALLEL TO 

45 ° FIBERS 


L I 

C) BETWEEN + 45/-45 PLIES 




d) in outside ( 45 °) ply 



Figure 37 (Continued). Crack/Delamination Growth in a Single-Edge 

Notched [± 45 / 0 ] Graphite/Epoxy Laminate at an 
Applied Axial llusile Stress of 105 MPa ( 15.2 ksi) 




a) BETWEEN 0/0 PLIES 


b) between 0/-45 PLIES 


Figure 38. Crack/Delamination Growth in a Single-Edge Notched [ ±45/0 ] g 
Graphite/Epoxy Laminate at an Applied Axial Tensile Stress 
of 300 MPa (44 ksi) . 
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c) BETWEEN -457+45 PLIES 


d) in outside (45°) ply 


Figure 38 (Continued). Crack/Delamination Growth in a Single-Edge 

Notched [±45/0] Graphite/Epoxy Laminate at 
an Applied Axial Tensile Stress of 300 MPa 
(44 ksi) . 
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a) BETWEEN 0/0 PLIES 


b) BETWEEN 0/-45 PLIES 


Figure 39. Crack/Delaminatino Growth in a Single-Edge Notched 

[±45/0] Graphite/Epoxy Laminate at an Applied Axial 
Tensile S Stress of 434 MPa (63 ksi) . 
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c) BETWEEN ~45/+45 PLIES 


d) in outside (45°) ply 


Figure 39 (Continued). Crack/Delamination Growth in a Single-Edge 

Notched [+45/0] Graphite/Epoxy Laminate at 
an Applied Axial Tensile Stress of 434 MPa 
(63 ksi) . 
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in predicting crack/ damage growth in composite laminates very similar to 
that which occurs in actual composites. 
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SECTION 8 


CONCLUSIONS AND SUGGESTIONS FOR FUTURE WORK 

A new fracture criterion has been developed for the fracture 
characterization of fiber-reinforced composites. The Griffith energy 
balance criterion has been modified to include energy absorption due to 
actual physical failure processes in the form of matrix yielding, matrix 
cracks, fiber breaks and fiber-matrix interface debonds, all of which 
are characteristic of the heterogeneous and anisotropic nature of these 
materials. The integrated micromechanical and macromechanical fracture 
criterion (IMMFC) proposed here is based on the assumption that micro- 
flaws in the form of broken fibers, matrix cracks and debonded fibers 
that exist in the composite for various reasons strongly influence the 
fracture process in these materials. Arguing that there is a relation 
between the size and density of the microflaws and the fracture response 
of the composite, parameters which are a measure of the size and density 
of the microflaws are introduced in the development of the IMMFC. 

A three-dimensional finite element micromechanics model has been 
presented to study the influence of these microflaws on the strength, 
stiffness and fracture toughness of a unidirectional lamina. The model 
has several advantages over two-dimensional transverse and longitudinal 
section models [58,59]. Even though the computational costs of using the 
three-dimensional micromechanics analysis are higher than for two- 
dimensional micromechanics analyses (the 3-D analysis uses 174 CPU 
seconds per solution increment on the CDC Cyber 760 computer, compared 
to 80 CPU seconds per solution when using the 2-D micromechanics analy- 
sis) , with the rapidly increasing capabilities of digital computers, the 



use of three-dimensional models such as in the present analysis appears 

to be the logical long term approach. 

In the present analysis, a local energy release rate in the 
presence of plasticity has been defined and used as a fracture criterion 
for the onset and growth of cracks. By including plasticity in the 
criterion, the analysis is applicable to most of the composite material 
systems currently in use. 

A crack growth simulation technique based on the virtual crack 
extension method has been developed for use in crack growth studies 
using a general three-dimensional finite element model. The accuracy of 
this crack growth technique has been established by applying it to 
predict the crack propagation in a centrally cracked aluminum plate. 

Several conclusions can be drawn from the results of the micro- 
mechanics analysis. For example, microcracks from broken fibers and 
fiber-matrix debonds grow in a self-similar fashion, whereas matrix 
cracks parallel to fibers are influenced by the presence of these fibers 
and change direction. This observation may be used to explain the 
experimentally observed zig-zag path of delamination cracks. The analy- 
sis can also be used to verify the double-cantilever beam test, in which 
a crack is introduced between unidirectional plies and propagated in the 
fiber direction. 

The macromechanics analysis, in which only delamination crack 
growth is allowed, is useful in studying the amount of delamination 

associated with different stacking sequences. 

In summary, the IMMFC and crack growth simulation techniques 
developed in the present analysis have been successfully used to predict 
the onset and growth of cracks/damage in the form of delaminations, 
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transverse cracks parallel to fibers, and through- the- thickness cracks 
in a single-edge notched [+45/0] graphite/epoxy laminate subjected to 
inplane tension normal to the notch. Even though the finite element 
grids used in the present analysis for both the micromechanics and the 
macromechanics analyses were shown to be adequate, a finer grid than the 
one used here will result in higher resolution. By considering different 
sizes and orientations of microflaws, a relation between microflaw 
parameters and the critical energy release rates eventually can be 
established. 

An extensive experimental study should also be initiated, to 
further establish the validity of the IMMFC. It should be possible to 
actually evaluate the material microflaw parameters, by observing and 
measuring the sizes of microflaws in actual composites. The critical 
energy release rates corresponding to the measured microflaw parameters 
could then be used for fracture characterization at the laminate level. 


Ill 
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APPENDIX A 


THREE-DIMENSIONAL, ELASTOPLASTIC, GENERALLY ORTHOTROPIC 
FINITE ELEMENT ANALYSIS 

The theory of the three-dimensional, elastoplastic, generally 
orthotropic finite element analysis developed by Monib and Adams [144] 
for the analysis of composite laminates is summarized here. 

The unidirectional composite lamina is assumed to be homogeneous 
and transversely isotropic. The three principal material axes are 
assumed as follows: the 1-axis, in the direction of the reinforcing 
fibers, the 2-axis, normal to the fibers but in the plane of the lamina, 
and the 3-axis, normal to the plane of the lamina (see Figure 2a of 
Section 3). The 2-3 plane is the plane of transverse isotropy. 

A quadratic yield condition similar to Hill’s yield criterion [152] 
is chosen in the form 

2f(o_) = F(o 2 - a 3 ) 2 + G(a 3 - c^) 2 + H(a 1 - a 2 ) 2 (A.l) 
+ 2L t 23 + 2M 4 + 2S r U ■ 1 

where F, G, H, L, M, and N are parameters characteristic of the current 
state of anisotropy. These parameters of anisotropy are allowed to vary 
with changes in temperature and/or moisture content. It can be shown 
[152] that 
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where a^, and a' are the tensile yield stresses in the 1, 2, and 3 

y y y 

directions of anisotropy. Also, if T ^ andx ^ are the yield 

stresses in shear with respect to the principal axes of anisotropy, then 
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(T 23> 2 
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(A. 3) 


The functional dependence of the parameters of anisotropy on temperature 
and moisture follows from Eqs. (A. 2) and (A. 3) when the yield stresses 
are expressed as functions of temperature and moisture content. 

The obvious association, implied by the term ’work-hardening, ’ 
between the work used to produce plastic flow and the hardening created, 
suggests the hypothesis that the degree of hardening is a function only 
of the total plastic work, and is otherwise independent of the strain 
path. The external work dW per unit volume done on the element during an 
infinitesimal increment of strain de^, with the continued loading of an 
element of the material is 


dW = a . ,de . . 
ij ij 


(A. 4) 



A part of this work 


dW = a..de®. (A. 5) 

e ij ij 


represents recoverable elastic energy; the remainder is the plastic work 
per unit volume, i.e.. 


dW = dW - dW = a. . (de . . - de?.) = a. .de?. 
p e xj 1 j ij ij ij 


(A. 6) 


P G 

where de.* = de.. - de. . is termed the plastic strain increment. The 
term dW^ is necessarily positive since plastic distortion is an irrever- 
sible process, and the plastic work is then 


W 

P 


* 


= a.. deK. 
ij ij 


(A. 7) 


In order for plastic work to be performed, the state of stress must 
be on the yield surface, i.e., the stress state must also satisfy the 
condition given by Eq. (A.l). To enforce this constraint, the Lagrange 
multiplier dA is used [153]. Then 

‘V E ij ' f( V dX1 ■ 0 <A - 8) 

which gives 


de? . = 9f dA 
X3 3c.. 


(A. 9) 


With the use of Eq. (A.l), a set of equations for the plastic strain 
increments can then be written as follows: 
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(A. 10) 

D * 

•?. = 2x..dX, 
13 13 

1 * 3 

(A. 11) 


where 

o* = [H (o x - 0 2 ) + G(o 1 - o 3 )]/(F + G + H) 
oj = [F(o 2 - o 3 ) + H(0 2 - 0 1 )]/(F + G + H) 

°3 = [G( °3 " + F( °3 ~ °2 )]/(F + g + h > 

T 23 = Lt 23 /(F + G + H) (A. 12) 

t*i = Mx 31 /(F + G + H) 
t* 2 = Nt 12 /(F + G + H) 

Separating the strains into elastic and plastic components gives 
de l = S ll dCT l + S 12 d °2 + S 13 d °3 + dG l 

de 2 = ^21 d °l + ^22 d °2 + ^23 da 3 + de 2 (A. 13) 

de 3 = S 31 da l + S 32 da 2 + S 33 da 3 + de 3 
and for the shear components 


dy 23 = S 44 dx 23 + dY 23 


dy 31 = S 55 dT 31 + dY 31 

(A. 14) 

dy 12 = S 66 dx 12 + dy 12 
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0 

where the S„ are the coefficients of the elastic compliance matrix [S ] 
in 

{de} » [S e ]{da> (A. 15) 

The inverse of Eq. (A. 15) is 

{da} = [C e ]{de> (A. 16) 

where [C 6 ] is the elastic stiffness matrix. If X , Y , etc., are the 

o o 

initial yield stresses, then according to the assumption of isotropic 
hardening above, X = hX Q , Y = hY Q , etc*, where h is a parameter increas- 
ing monotonically from unity that expresses the amount of hardening. The 

anisotropic parameters must then decrease in accordance with Eq. (A. 2) 
2 

as F = F Q /h , etc. The manner in which h varies with strain can be 
explained by analogy with the isotropic theory due to von Mises. Let 
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+ F + G + H 

be a nondimensional measure of the equivalent stress a. By analogy with 
the von Mises criterion for isotropic materials. Hill [152] suggested 
that if there is a functional relation between o and the work W (this is 
yet to be demonstrated by experiment) , there must be one between a 
and the effective (or equivalent) strain, / de, as defined by 
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where de„ and dy.^ are given by Eqs. (A. 13) and (A. 14). This is the 
analogue of the equivalent stress-equivalent strain curve for isotropic 
materials, the area under which is equal to the work per unit volume. 
Accordingly, 


dW = o(de - de 6 ) = ode* 3 


(A. 19) 


But, from Eq. (A. 5) 


dW = a, de, + . - 
P 11 


,. + x 23 dy 23 + 


(A. 20) 


Substituting for de?.. and dy?^. from Eq. (A. 10) into Eq. (A. 20) yields 


2 -2 

dW = -r- a dX 
P 3 


(A. 21) 


If an effective stress-effective plastic strain curve is then 
constructed, the slope of such a curve at any point will be 


H’ = 


do 

diP 


{A. 22) 


from which 


j-P do 
de = W 


(A. 23) 
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Now, substituting for de P in Eq. (A. 19), and equating the result to Eq. 

(A. 21) since both are equal to the plastic work per unit volume dW , 

P 


do = 4 5 dX = dW 
H J p 


(A. 24) 


Rearranging 


2 _ _ 4 _ 

-j odo ~ ~g a H'dX 


(A. 25) 


The left-hand side of Eq. (A. 25) is recognized as the differential of 
Eq. (A. 17), defining a. Thus, 


2 -j- 

3 ado = ? 


T g"+"h < [ [ -G ( a 3 “ a 1 ) da 1 + HC 0 ! - a 2 ^ da i^ 


+[F(a 2 - 03)da 2 - H(a 1 - o 2 )da 2 ] 
+[-F(a 2 - 03)do3 - 0(03 - 03^03] 


+ 2Lx 23 dT 2 3 + 2MT 31 dT 31 + T2NT 12 dT 12 


(A. 26) 


With the use of the definitions of Eq. (A. 12), Eq. (A. 26) can be re- 
written as 

2 _ _ * * * 

-3 odo = OjdOj^ + o 2 da 2 + °3 da 3 

+ 2x 23 dT 2 3 + 2x3 1 dx 31 + 2x 12 dx 12 (A. 27) 

To arrive at a relation between stress and strain, Eqs. (A. 13) and 
(A. 14) must be inverted. Rewriting these equations in matrix form, 
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{de} = [S e ]{d a} + {de P } 
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(A. 28) 



C ll°l + C 12°2 + C 13 a 3 


C 12 CT 1 + C 22 a 2 + C 23 a 3 


C 13°l + C 23°2 + C 33°3 


{A} = < 


2C 44 T 23 


2C 55 T 31 


2C 66 t 12 


(A. 35) 


Equating Eqs. (A. 25) and (A. 27) yields 


j o 2 H'dX = 


k k k 

°l^ a l + a 2^ a 2 + a 3^ a 3 


+ 2x 23 dr 23 + 2x 31 dT 31 + 2x 12 dx 12 


(A. 36) 


Now, substituting for da in Eq. (A. 36) from Eq. (A. 34), and solving for 
dX 


dX 


A^de^ 4* A 2 de 2 i A 3 ds 3 i A^dy 23 4* A^dy 3 ^ A^dy^ 2 


(A. 37) 
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where A^(i =1, ...» 6) are elements of {A}, and 
B = 3 2 H' + AjO* + A 2 o* + A 3 o* 

+ 2A^t 33 + 2A^t 33 + 2A^t^2 (A. 38) 

Finally, substituting for dX from Eq. (A. 37) into Eq. (A. 38) yields the 
desired form for the stress-strain relation, viz, 

(do) = [ C P ] { de } (A. 39) 

where 



is the plastic stiffness matrix. 
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For the 20-node quadratic isoparametric brick element, the element local 
coordinates and node numbers used in the present formulation are as 
shown in Figure A.l. 





Figure A.l. Node Point Numbering System for the 20-Node Quadratic 
Isoparametric Brick Element. 


The shape functions are: 
For the corner nodes 


“i ■ 5 (1 + £ o> 


t-L 


i = 1-8 


(A. 41) 
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For the midside nodes 


N. = i (1 - C 2 ) a + n o ) (1 + c 0 ) i = 9-12 

i (1 - n 2 ) (1 + 5 q ) (1 + 4 Q ) i = 13-16 (A. 42) 

N. = ^ (1 - 4 2 ) (1 + ? Q ) (1 + n o ) i = 17-20 

where 

£ o = ec i ’ n o = nn i ’ and C o = «i 

The relation between derivatives in the local and global coordinate 
systems is given by the chain rule of differentiation as 



where commas denote partial differentiation and [J] is the tranformation 
Jacobian matrix, which can be given as 


‘ BNj/H 3N 20 /35 I y 1 

[J] = 3N x /3n 3N 20 /3n x 2 y 2 z 2 (A. 44) 

. s V 8t 3K 2 o / s? J : : : 

• • • 

x 20 y 20 Z 20 
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The strains at a point within an element can be written in terms of 
element nodal displacements as 


{e} = [B] {6 } 

where { 5}^ = [u, v, w, v 2 ’ w 2 *** u 20’ V 20* W 20^ 
vector. 

The element stiffness matrix is given by 


[K] = f [B] T [C][B] d(Vol) 

Vol 

The C coefficients in the generally orthotropic 
Eq.(A.46) are listed below: 


C = \J 1 + U 2 cos (20) + U 3 cos (40) 
C 2 = U 4 + U 3 cos (40) 


2 2 

C^ 3 *» C^cos 0 + C2 3 sin G 


'14 


'16 


'22 


U 


23 


'24 


'26 


c i5 - 0 

U sin (20) - U 3 sin (40) 

cos(20) + U 3 cos (40) 

2 2 
C^sin 0+ C23COS G 

C 25 - 0 

\ U 2 Sin 20 + U 3 sin(40) 


'33 


= C 


34 

^36 = 


'44 


33 


C 35 = 0 

= Cj^sin© cos0 - C 22 Sin 0 cos 0 

2 2 
= C..cos 0+ C c( .sin 0 
44 55 


= C^sin© cos© + C^sin© cos© 


(A. 45) 

is the displacement 
(A. 46) 

stiffness matrix in 


(A. 47) 
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E 46 - 0 

2 2 

C ct . = C,,sin 0+ C cc cos G 
55 44 55 

E 56 ‘ 0 

E 66 * U 5 ’ U 3 cos (48) 

where 

U x = 1/8 (3 C 1:l + 3C 12 + 4C 66 ) 


U 2 " 

1/2 

< c n - 

c 22 > 



U 3 = 

1/8 

< c n + 

C 22 

- 2C 12 

- «66 

U 5 = 

1/8 

< c u + 

C 22 

' 2C 12 

+ 4C 66 

C 11 

= (1 

- v 23 

V 22 ) 

VE 11 


C 22 

= (1 

" V 31 

V 13 ) 

ve 22 


C 33 

= (1 

_ V 12 

V 21 ) 

VE 33 



C 12 = (V 21 + V 23 V 31 } VE 11 = (V 12 + V 13 V 32 ) VE 22 

C 13 = (V 31 + V 21 V 32 } VE 11 = (V 13 + V 23 V 12 } VE 33 

C 23 = (V 32 + V 12 V 31^ VE 22 = (v 23 + V 21 V 13^ VE 33 

C 44 = G 23 
C 55 = G 31 
C 66 = G 12 

and 

v = d - v 12 v 21 - v 23 v 32 _ v 13 v 31 _ v 12 V 23 V 31 ) 

\ 

The Richard-Blacklock three-parameter model [148] has been 
represent the nonlinear stress-strain response of the material, 
model is in the form: 

a - Ee/ [1 + |Ee/o 0 | n ] 1/n 


(A. 48) 


(A. 49) 

used to 
The 

(A. 50) 
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where o and n are two independent parameters, and E is the initial 

o 

slope of the stress-strain curve. Since the shape of an effective 
stress-effective strain curve is similar to a tensile stress-strain 
curve, a similar equation for the effective stress-effective strain can 
be written as 

o = E e/[l + |E e/3 o | n ] 1/n (A. 51) 

where a is the effective stress and Z is the effective strain, as 
defined by Eqs. (A. 17) and (A. 18), respectively. The two independent 
parameters o q and n, together with the third parameter E, which is the 
initial slope of the curve, are selected to best fit the experimental 
data. 

The tangent modulus H' of the effective stress-effective plastic 
strain curve is related to E as [154] 


H' 



(A. 52) 


where E^ is found by differentiating Eq. (A. 51) with respect to £. The 
resulting equation is 

1+n 

E t = E/[l + |E £/aJ n ] n (A. 53) 

Since the material properties depend on temperature in the case of 
metal matrix composites and on both temperature and moisture in the case 
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of polymetric matrix composites, the material properties are represented 
by a second order polynomial of the form 

P(T,M) = C l T 2 + C 2 M 2 + C 3 TM + C 4 T + C 5 M + C & (A. 54) 

where P is the functional material property of interest, T is temper- 
ature, M is moisture content, and through are regression coef- 
fic Lents for that property. The parameters in Eq. (A. 51), viz E, a Q and 

n, which can also vary with temperature and moisture, are thus also 
expressed in the same polynomial form of Eq. (A. 54). 
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APPENDIX B 


Crack Growth Simulation Technique 

The crack growth simulation in a three-dimensional elastoplastic 
finite element analysis, as developed in Section 4, is further explained 
here. Implementation of this crack growth simulation technique in the 
actual finite element computer code is also discussed. 

The technique uses the critical strain energy release rate 
(fracture toughness) of the material as the criterion for the onset and 
growth of a crack. That is, at any stage of loading, if the computed 
strain energy release rate for a virtual crack extension Aa exceeds the 
fracture toughness of the material, the original crack is assumed to 
have extended by an amount Aa. 

The strain energy release rate in the presence of plasticity, I, is 
computed using the rate of change of compliance method. It must be noted 
that the symbol I is used to denote strain energy release rate in the 
present analysis, to distinguish it from the symbol G which is generally 
used to represent the strain energy release rate in linear elastic 
fracture mechanics. 

Making use of the observations made in Section 4 relative to the 
stiffness derivative method for both fixed load and fixed grip 
conditions, that any small extension of the crack in a finite element 
model affects only the relatively few elements near the crack tip, the 
strain energy release rate in the present analysis is computed in terms 
of local compliance changes. 

The organization of the three-dimensional, elastoplastic, generally 
orthotropic finite element analysis computer code (WY03D) with crack 



propagation capability, as developed for use in the present analysis, is 
shown in Figure B.l. 



Figure B.l. Block Diagram of the WY03D Computer Code 
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The bulk of the input data regarding the model geometry, material 
properties, and loading conditions are read into the INPUT subroutine 
following the general format given in Reference [144]. For crack 
propagation, the following additional variables should be read as input: 

NPCR: Crack Propagation Flag. It must be equal to 1 for crack 

propagation analysis. 

NNCT: Number of Nodes near a Crack Tip. 

NNCTN : Node Number of (brack Tip Nodes. 

NELCN : Elements surrounding a crack tip node as shown in 
Figure B.2. 



Crack Tip Node 


Figure B.2. Numbering Order of Elements Near a Crack Tip Node 
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The step by step procedure adapted for crack growth simulation in 

the computer code is as follows: 

Step 1: During the first load increment, the total stiffness 

coefficients (used for computing the strain energy release 
rates) at all crack tip node points in all nine modes of crack 
extension (refer to Figures 5, 6, and 7) are evaluated in 
subroutines STRESS, CRACK, and STIFF, as explained in Section 
4 using a two-dimensional example. The stiffness coefficients 
K are then stored in an array for later use. 

Step 2: The three incremental force components F x , F^, and F^ at all 

crack tip nodes are also stored in an array (this is done in 
subroutine CRACK1) . In all of the subsequent load increments, 
the incremental force components F^, F^, F z at the crack tip 
nodes are computed and added to the previous incremental 
values to obtain total force components (in CRACK1) . 

Step 3: In subsequent load increments, the crack is virtually extended 

to the next nearest grid points in all possible directions (in 
both opening and shearing modes) and the reduced stiffness 
coefficients K are evaluated as explained in the two- 
dimensional example in Section 4. It must be noted that there 
may not be two values of K for each mode, depending on which 
set of elements were used to compute K (e.g.. Elements 1 and 
2, or 3 and 4 in Figure 9). 
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Step 4: 


The strain energy release rates I are calculated in all nine 


inodes of crack extension using the expressions given in Eq. 
(29). (Equation 29 gives I values for Modes M T , M , and 

M ZIII* F ° r Modes Mxi’ M XIII* ^1’ ^II* and ^III* 

similar expressions can be derived by cylically changing the 

subscripts.) In an actual computation, there are two values of 

I for each mode, corresponding to two values of the reduced 

stiffness coefficients. The higher of the two computed strain 

energy release rates is taken as the strain energy release 

rate in that particular mode. 

Step 5: In subroutine CRACK2 the computed values of the strain energy 

release rates are compared with the critical energy release 
rates (the fracture toughness values for the nine modes of 
fracture, which are read as material properties in subroutine 
INPUT). If the computed strain energy release rate is higher 
than the critical strain energy release rate, the crack tip 
node point is split into two node points. A new node number 
is then assigned to one of the nodes and the element 
connectivity matrix is updated to include new nodes. Reaction 
forces are than applied to the old and the new node points. 
The crack propagation flag (NPCR) is set equal to 2, to check 
for any further crack growth due to reaction forces in the 
same load increment. If there are no further crack extensions, 
the NPCR is set back to 1 and loading is resumed. 
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